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ABSTRACT 


Parameter  estimation  methods  were  used  to  obtain  estimates  of  stability  and  control  derivatives 
for  the  Marine  Corps  BQM-147  Unmanned  Air  Vehicle.  The  results  from  a  simple,  PC-based  linear 
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was  used  to  assess  the  accuracy  of  the  estimates  for  both  methods.  The  bounds  were  high  for  both 
the  longitudinal  case  and  the  lateral-directional  case  due  to  the  limited  maneuvers  tested,  high  levels 
of  noise  in  the  same  general  frequency  range  as  the  control  input,  and  the  lack  of  body-angle  data. 
The  linear  model  failed  to  provide  estimates  for  the  lateral-directional  case.  Though  the  results  may 
be  used  as  starting  points  for  a  dynamic  model  of  the  aircraft  it  is  recommended  that  the  flight  test 
procedures  be  modified  to  address  the  issues  raised  concerning  noise,  recorded  signals  and  the  need 
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I .  INTRODUCTION 

Symmetrical-wing,  delta-planform  (SYMDEL)  unmanned  air 
vehicles  (UAVs)  have  been  in  development  and  have  flown  in  a 
variety  of  mission  profiles  since  the  early  1970 's.  This 
delta-planform,  symmetrical-airfoil  configuration  was  designed 
to  provide  equal  performance  in  range,  speed  and  stability  in 
either  an  upright  or  inverted  attitude  [Ref.  l:p.  23].  With 
this  attribute  the  vehicle  can  readily  accommodate  a  variety 
of  payloads  which  makes  it  ideally  suited  for  specific 
military  applications  or  as  a  testbed  vehicle  for  new 
component  technologies. 

The  EXORONE  UAV  was  developed  in  the  early  1980 's  at  the 
Johns  Hopkins  University,  Applied  Physics  Laboratory  (APL) . 
The  vehicle  design  is  a  descendant  of  the  earlier  SYMDEL  4 
with  a  constant-chord  section  inserted  into  the  middle  of  the 
delta  planform  to  provide  an  accessible  avionics  bay  [Ref. 
l:p.  25].  The  original  APL  version  of  the  EXDRONE,  referred 
to  as  the  baseline  configuration  throughout  this  document,  was 
intended  to  meet  a  United  States  Marine  Corps  (USMC) 
requirement  for  a  UAV  to  carry  an  expendable  VHP 
communications  jammer.  Takeoff,  climb  out  and  turn  to  the 
desired  heading  are  accomplished  by  radio  control.  Once  the 
mission  heading  is  achieved,  a  radio  control  command  switches 
control  of  the  vehicle  to  the  autopilot.  The  autopilot  then 
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controls  the  vehicle  for  the  duration  of  the  mission.  The 
name  EXDRONE  implies  that  the  vehicle  was  intended  to  be 
operated  autonomously  and  to  be  expendable  at  the  end  of  its 
mission;  however,  radio  control  can  be  reacquired  if  necessary 
and  the  vehicle  can  be  recovered. 

In  preliminary  flight  tests  the  baseline  configuration  of 
the  EXDRONE  was  found  to  be  difficult  to  handle  for  less 
experienced  pilots.  The  vehicle  had  weak  stability  and 
control  characteristics  in  low-speed  flight  and  had  a  tendency 
to  depart  near  stall  [Ref.  2:p.  1].  The  primary  mission  of 
the  EXDRONE  is  to  dash  outbound  at  a  speed  of  about  100  miles 
per  hour  to  a  predesignated  target  area  and  then  to  loiter  in 
the  area  and  jam  enemy  VHP  communications.  With  low  speed 
flight  being  an  equally  important  part  of  the  EXDRONE  mission, 
the  USMC  UAV  Project  Office  contracted  with  the  NASA  Langley 
Research  Center,  Flight  Dynamics  Branch,  to  conduct 
exploratory  wind  tunnel  and  flight  tests  in  order  to  determine 
the  aerodynamic  characteristics  of  the  vehicle  and  to  provide 
modifications  for  aerodynamic  improvements  [Ref.  2:p.  1].  As 
a  result  of  the  wind  tunnel  study  several  modifications  were 
made  to  the  baseline  configuration  in  order  to  increase  pitch 
control  for  nose-up  trim,  improve  the  stall  departure 
characteristics,  minimize  the  adverse  yaw  characteristic  and 
reduce  the  minimum  level  flight  speed  as  a  means  of  increasing 
endurance.  It  was  this  modified  configuration  EXDRONE  that 
became  the  subject  of  the  flight  test  portion  of  the  NASA 
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Langley  study.  A  detailed  description  of  the  wind  tunnel 
tests  on  both  the  baseline  and  modified  configurations  can  be 
found  in  Reference  2. 

The  modifications  to  the  wing  tips  and  control  surfaces 
that  are  described  in  Reference  2  were  included  in  the 
production  contract  for  the  BQM-147A.  The  following 
additional  modifications  are  also  a  part  of  the  production 
conf iguration . 

1)  Position  and  orientation  of  the  engine 

2)  Shape  of  the  nose/cowl 

3)  Number  and  type  of  skids 

In  October  of  1991,  Developmental  Testing  (DT)  was 
conducted  on  the  BQM-147A  at  the  Dugway  Proving  Grounds,  Utah. 
The  purpose  of  the  tests  was  to  determine  the  suitability  of 
the  vehicle  in  the  VHF  communications  jamming  mission  and  to 
determine  the  effects  of  the  communications  jammer  payload  on 
the  handling  qualities  of  the  airframe.  The  BQM-147A 
control  surfaces  are  manipulated  by  electrical  servos  that  are 
controlled  by  an  autopilot.  The  autopilot  receives  command 
signals  from  either  a  modular  radio  control  unit  for  direct 
pilot  input  or  from  an  onboard  microprocessor  containing  a 
pre-programmed  mission  scenario  for  the  autonomous  mode.  The 
autopilot  processes  the  error  signal  from  a  rate  gyro  to 
augment  roll  damping  in  flight  [Ref.  3:p.  20].  During  DT  the 
autopilot  had  to  be  fine  tuned  by  hand  in  a  hit  or  miss 
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fashion  because  the  vehicle  control  laws  had  not  been  well 
defined  during  the  pre-production  work.  It  was  at  this  point 
that  the  USMC  UAV  Project  Office  decided  to  fund  a  series  of 
flight  tests  from  which  time  history  data  could  be  used  to 
determine  the  aerodynamic  characteristics  of  the  aircraft. 

Research  being  done  at  the  Naval  Postgraduate  School 
(NPS) ,  Department  of  Aeronautics  and  Astronautics,  using  UAV 
testbeds  has  been  well  received  by  the  UAV  Joint  Project 
Office.  Specifically,  aerodynamic  parameters  obtained  by 
numerical  methods  and  by  wind  tunnel  tests  [Refs  4  and  5]  have 
been  used  by  the  Target  Simulation  Lab,  Naval  Air  Warfare 
Center,  Weapons  Division  (formerly  Pacific  Missile  Test 
Center)  ,  Pt.  Mugu,  CA.,  in  the  development  of  a  math  model  for 
the  Pioneer  UAV  flight  training  simulator.  More  recent  thesis 
work  by  U.S.  Navy  LCDR  Robert  Graham  incorporated  a  personal 
computer  (PC)  based  parameter  estimation  capability  into  the 
NPS  research  program  [Ref.  6].  Being  advised  of  this  desktop 
parameter  estimation  capability,  the  USMC  UAV  Project  Office 
expressed  an  interest  in  an  application  of  this  routine  using 
measured  data  obtained  from  the  next  series  of  flight  tests. 

The  objectives  of  the  application  listed  in  Reference  7  are: 

1)  To  document  the  performance  of  the  aircraft. 

2)  To  reprogram  the  autopilot  on  the  production  aircraft. 

3)  To  form  a  basis  for  suggesting  changes  to  the  production 
aircraft  which  will  improve  stability  and  performance. 

4)  To  develop  a  PC-based  flight  simulator  which  will  be 
used  to  train  operators. 
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Aircraft  parameter  estimation  routines  are  used  to 
estimate  aircraft  stability  and  control  derivatives  from 
actual  time  history  data.  The  stability  and  control 
parameters  extracted  from  these  routines  provide  an  accurate 
model  of  the  actual  vehicle  and  can  be  used  in  analyzing  the 
effects  of  proposed  control  system  changes  on  aircraft 
stability  and  handling  qualities.  Parameters  obtained  from 
observed  data  can  also  be  used  to  provide  an  accurate  model  of 
the  vehicle  for  greater  fidelity  in  flight  simulators. 

The  purpose  of  this  study  was  to  use  the  flight  test  time 
history  data  found  in  Appendix  B  in  a  parameter  estimation 
model  to  estimate  both  the  longitudinal  and  lateral- 
directional  stability  and  control  derivatives.  Two  parameter 
estimation  models  were  used.  The  first  was  an  application  of 
a  PC-based  parameter  estimation  program  by  LCDR  Robert  Graham 
[Ref.  6].  This  model  was  selected  because  it  is  readily 
available,  it  operates  in  a  familiar  DOS  environment  and  it 
generally  provides  good  results  when  compared  with  more 
complex  mainframe  programs  [Ref.  8].  The  disadvantage  of  this 
first  parameter  estimation  routine  is  that  it  is  a  linear, 
maximum  likelihood  estimator  which  assumes  that  the  actual 
system  satisfies  the  assumptions  that  were  made  in  the 
linearization  of  the  model.  The  second  model,  named  pEst,  is 
an  interactive  program  for  parameter  estimation  in  non-linear 
dynamic  systems.  Incorporation  of  the  pEst  routine  into  the 
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ongoing  flight  test  research  at  NPS  provides  a  capability  to 
model  the  dynamic  behavior  of  aircraft  that  cannot  be 
appropriately  modeled  using  the  assumptions  inherent  in  the 
shorter,  PC-based  MMLE-3  code.  Because  of  suspected  non- 
linearities  in  the  more  extreme  areas  of  the  flight  envelope 
of  the  EXDRONE,  the  linear  model  may  not  be  sufficient  and  the 
pEst  routine  might  better  capture  the  behavior  of  the  system. 
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II.  BACKGROUND 


A.  HISTORICAL  REVIEW 

The  concept  of  using  remotely  piloted  vehicles  (RPV's)  and 
unmanned  aerial  vehicles  (UAV's)  in  a  tactical  combat  role  had 
limited  acceptance  within  the  military  community  from  the  time 
of  its  first  proposal  in  1915  until  early  in  the  last  decade. 
The  RPV  is  controlled  by  a  human  operator  by  means  of  a  radio¬ 
control  guidance  link  using  either  direct  operator  visual 
contact  or  onboard  video  imagery  to  sense  flight  parameters. 
Initial  RPV  testing  was  limited  by  pre-WW  II  radio-control 
technology  but  improvements  in  guidance  technology  led  to 
successful  developmental  flight  tests  of  ordnance  carrying 
RPV’s  prior  to  the  attack  on  Pearl  Harbor  in  December  1941. 
There  was  limited  use  of  these  lethal  RPV's  throughout  the 
latter  stages  of  World  War  II  and  the  Korean  Conflict. 
Widespread  use  of  these  platforms  and  full-scale  development 
of  follow-on  vehicles  with  greater  combat  potential  were  met 
with  opposition  by  military  planners  who  preferred  to  use 
scarce  resources  on  more  expensive  and  ostensibly  more 
technically-challenging  guided-raissile  programs  [Ref.  10]. 

Advancements  in  radio-control  and  digital  communication 
technologies  coupled  with  an  operational  need  to  defeat  more 
modern  military  anti-air  warfare  weapons  led  to  an  expanded 
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role  for  a  new  generation  of  non-lethal  UAV's.  The  UAV 
differs  from  the  RPV  in  that  the  guidance  signal,  from  either 
a  human  operator  via  radio-control  uplink  or  from  on-board 
navigation  subsystems,  is  input  to  an  auto-pilot  which 
generates  the  appropriate  control  deflection.  UAV's  are 
generally  known  for  their  ability  to  operate  autonomously  over 
the  entire  range  of  a  preprogrammed  mission  profile.  During 
the  Vietnam  conflict  the  Ryan  Model  147  Firebee  was  used  for 
photo  reconnaissance,  electronic  intelligence  gathering, 
communications  intelligence  and  covert  psychological  warfare 
missions  [Ref.  ll:p.  20].  In  spite  of  the  tremendous  success 
of  UAV's  in  Southeast  Asia,  continued  support  for  further  UAV 
development  waned  in  the  wake  of  higher  priority  military  and 
domestic  programs.  Not  until  the  unprecedented  and  highly 
successful  use  of  UAV's  by  the  Israelis  in  the  Bekaa  Valley  in 
1982  did  U.S.  military  planners  become  convinced  of  the 
tactical  usefulness  on  non-lethal  UAV's  [Ref.  12]. 

Recognizing  the  need  for  common  and  interoperable  systems. 
Congress  in  1988  directed  the  Department  of  Defense  (DoD)  to 
consolidate  the  management  of  non-lethal  UAV  programs.  The 
DoD  established  the  UAV  Joint  Project  Office  (UAV  JPO)  and 
appointed  the  United  States  Navy  as  the  executive  service. 
The  mission  of  the  Joint  Program  Office  is  to  expeditiously 
field  quality  UAV  systems  which  provide  a  significant  tactical 
advantage  to  operational  commanders.  The  objective  of  the  UAV 
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JPO  is  to  field  a  family  of  interoperable  systems  that 
optimize  commonality  [Ref.  13:p.  7]. 

During  Operation  Desert  Storm,  UAV's  provided  tactical 
commanders  an  operational  capability  that  was  not  previously 
available.  The  Marines  operated  three  types  of  UAV's  in 
various  mission  roles  that  included  real  time  reconnaissance, 
surveillance,  target  acguisition,  spotting  and  bomb  damage 
assessment.  When  operation  of  the  Pioneer  UAV  was  jeopardized 
by  parts  shortages  and  while  the  more  fragile  Pointer  was 
unable  to  fly  in  winds  greater  than  15  knots,  the  Marine  Corps 
augmented  their  UAV  operations  by  reassigning  a  number  of 
EXDRONE 's  from  a  research  and  development  program  to  be 
operationally  deployed  from  airfields  near  the  Saudi  Arabia- 
Kuwait  border  [Ref.  14].  These  vehicles  were  reequipped  with 
color  video  cameras  and  were  used  to  find  attack  routes 
through  enemy  defenses. 

This  demonstrated  capability  of  UAV's  to  perform  multiple 
functions  and  the  ability  to  provide  near  real-time 
intelligence  in  a  high-threat  environment  without  loss  of  life 
has  affirmed  the  commitment  by  the  Services  to  integrate  UAV's 
into  the  force  structure.  The  1992  DoD  Unmanned  Aerial 
Vehicles  Master  Plan  [Ref.  13]  indicates  that  the  planned 
procurement  of  UAV's  is  in  excess  of  2300  vehicles.  The  core 
strategy  of  the  UAV  JPO  is  the  establishment  of  a  family  of 
UAV  systems  which  indicates  that  these  systems  must  allow  for 
growth  in  performance  and  must  readily  accommodate  changing 
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payloads.  As  a  part  of  this  core  strategy  the  JPO  has 
initiated  a  Very  Low  Cost  (VLC)  UAV  program.  These  VLC 
systems  are  used  to  demonstrate  and  evaluate  the  utility  of 
UAV's  in  tactical  units. 

In  November  1991,  the  UAV  JPO  awarded  a  contract  to  BAI 
Aerosystems,  Inc.  for  the  production  of  110  EXDRONES.  At  the 
present  time,  operational  test  and  production  for  the  Marine 
Corps  Jammer  program  is  not  funded  and  the  EXDRONE  vehicles 
are  being  used  in  the  JPO  VLC  program  [Ref.  13:p.  49].  If  the 
production  option  in  the  BAI  contract  is  exercised,  it  is 
expected  that  as  a  result  of  the  VLC  program  a  new  effort  to 
improve  range  and  payload  capacity  will  begin.  Having 
accurate  estimates  of  stability  and  control  derivatives  will 
allow  designers  to  incorporate  good  flying  qualities  from  the 
earliest  stages  of  design  and  avoid  repeating  errors  from 
earlier  versions  of  the  model. 


B.  PARAMETER  ESTIMATION 

In  Reference  15  Ljung  points  out  that  the  construction  of 
a  model  from  observed  data  involves  three  basic  entities; 

1)  The  data  record. 

2)  The  set  of  candidate  models. 

3)  A  rule  by  which  candidate  models  can  be  assessed  using 
the  data. 
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In  an  application  to  aircraft  stability  and  control,  the 
flight  test  maneuvers  should  be  designed  so  that  the  data 
collected  will  provide  the  maximum  value  to  the  analyst. 
Conseguently  it  essential  that  the  flight  test  program  have  a 
clear  set  of  objectives  so  that  the  flight  test  maneuver  will 
cover  the  area  of  the  envelope  that  is  to  be  explored. 
Instrumentation  must  be  considered  so  that  the  appropriate 
parameters  are  measured  to  a  degree  of  accuracy  that  will  be 
suitable  to  the  objectives  of  the  program.  The  set  of 
candidate  models  must  also  be  considered.  The  equations  of 
motion  derived  from  analysis  of  forces  and  moments  acting  on 
a  typical  aircraft  are  well  defined.  For  most  applications  a 
linearized  set  of  the  equations  is  adequate  for  analyzing 
aircraft  stability  and  control.  If  a  linear  model  is  all  that 
is  available  for  the  analysis,  then  only  small-  perturbation 
maneuvers  should  be  included  in  the  flight  test.  If  envelope 
expansion  is  the  objective  of  the  flight  test,  then  maneuvers 
that  gradually  expand  the  envelope  may  provide  data  that 
indicate  a  non-linear  model  is  necessary.  The  rule  by  which 
the  models  are  assessed  is  generally  determined  by  how  well 
the  model  is  able  to  reproduce  the  measured  data.  For 
example,  the  modified  maximum  likelihood  program  uses  an 
output-error  approach  that  minimizes  the  difference  between 
the  measured  parameters  and  the  model  estimate  of  the 
parameters  in  successive  iterations.  This  construction  of  a 
mathematical  model  of  a  dynamic  system  from  observed  data  is 
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know  as  System  Identification  [Ref.  15:p.  1]. 

Estimation  of  aircraft  stability  and  control  derivatives 
from  flight  test  time  history  data  is  one  of  many  practical 
applications  of  the  basic  system  identification  techniques. 
Of  particular  note  in  the  aeronautics  industry  is  the  work  of 
Lawrence  Taylor,  Kenneth  Iliff ,  Richard  Maine  and  James  Murray 
of  the  NASA  Ames  Research  Center,  Dryden  Flight  Research 
Facility.  Their  approach  to  the  estimation  of  aircraft 

stability  and  control  derivatives  through  the  application  of 
system  identification  techniques  can  be  found  in  References  9, 
17  and  18.  The  method  of  application  can  be  summer ized  in  the 
following  steps: 

1)  Formulate  the  dynamic  equations  that  describe  the  model. 

2)  Identify  from  the  form  of  the  dimensional  derivatives 
the  parameters  whose  values  are  unknown. 

3)  Collect  flight  test  data. 

4)  Infer  the  values  of  the  unknown  parameters  by  adjusting 
the  values  iteratively  until  the  computed  response  best 
matches  the  measured  response. 

Throughout  the  1980 's  the  most  widely  accepted  aircraft 
parameter  estimation  routine  was  the  Modified  Maximum 
Likelihood,  version  3,  written  by  Iliff  and  Maine.  The  MMLE-3 
program  evolved  from  an  earlier  Newton-Raphson  parameter 
estimation  routine  by  Taylor  and  Iliff  when  a  more  versatile 
program  capable  of  handling  larger  amounts  of  data  was 
required  by  the  industry.  When  applying  system  identification 
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techniques  to  a  real  world  system,  modeling  errors  are 
unavoidable  since  it  is  likely  that  the  real  system  does  not 
completely  satisfy  all  of  the  assumptions  on  which  the  model 
was  developed.  The  existence  of  unmeasured  inputs  along  with 
these  modeling  imperfections  perturb  the  states  and  to  at 
least  some  degree  invalidate  the  results  of  the  identification 
algorithms  [Ref.  15].  The  advantage  of  MMLE-3  is  that  the 
algorithm  is  capable  of  applying  the  system  identification 
techniques  in  the  presence  of  state  noise  and  thus  produces  a 
more  optimistic  error  bound  estimate.  Full  details  of  the 
MMLE-3  program  can  be  found  in  References  17  and  19.  The 
minimization  routines  and  the  criterion  used  in  a  personal 
computer  based  adaptation  of  MMLE-3  are  discussed  in  more 
detail  in  later  sections. 

For  many  aircraft  applications  the  linearized  equations  of 
motion  provide  a  reasonably  accurate  model  of  the  system  for 
small  disturbances.  For  some  applications  however,  the  linear 
model  is  not  sufficient  and  a  more  accurate  non-linear  model 
is  required.  Flight  testing  at  the  Dryden  Flight  Research 
Facility  highlighted  the  limitations  of  the  linear  model  when 
applied  to  maneuvers  at  extreme  flight  conditions  and  with 
some  unique  aircraft  configurations  that  exhibited  non-linear 
dynamic  behavior.  In  1987  Murray  and  Maine  introduced  a 
program  for  parameter  estimation  in  non-linear  dynamic 
systems.  The  capabilities  of  the  pEst  program  are  described 
in  Reference  9. 
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Historically  the  wind  tunnel  has  served  as  the  primary 
source  of  aeronautical  data  for  flight  vehicles.  The 
advantage  of  the  wind  tunnel  is  that  flow  conditions  can  be 
accurately  controlled,  but  other  influences  such  as  wall 
interference  are  difficult  to  account  for.  The  wind  tunnel  is 
also  limited  in  its  ability  to  achieve  dynamic  similarity  with 
actual  flight  test  conditions.  There  are  numerous  reasons  why 
wind  tunnel  data  alone  are  not  sufficient  and  often  the 
results  are  validated  by  flight  test.  In  Reference  18,  Maine 
and  Iliff  point  out  that  there  is  a  tendency  in  practice  to 
emphasize  positive  results  and  to  down  play  the  role  of  flight 
test  in  validating  wind  tunnel  data  or  analytical  predictions. 
A  discussion  of  the  errors  that  are  found  and  corrected 
because  of  the  flight  test  are  often  omitted  from  the 
published  reports.  Maine  and  Iliff  go  on  to  say  that  the 
utility  of  stability  and  control  parameter  estimation  can  best 
be  evaluated  by  comparing  the  predicted  data  prior  to  the 
first  flight  test  with  the  best  estimates  from  combined  flight 
data  and  predictions  at  the  conclusion  of  the  flight  test 
program  [Ref.  18 :p.  4].  The  task  then  is  to  take  maximum 
advantage  of  previous  information,  both  analytical  and 
experimental,  and  to  use  the  estimation  routines  to  fill  the 
gaps.  The  parameter  estimation  routines  used  in  this  study 
use  a-priori  knowledge  of  the  system  under  investigation  to 
provide  to  the  user  the  most  accurate  estimates  of  aircraft 
stability  and  control  derivatives. 
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III. 


DEVELOPMENT  OF  THE  MATHEMATICAL  MODEL 


The  analysis  of  a  dynamic  system  begins  with  the 
derivation  of  a  mathematical  model  that  is  reasonably 
accurate.  Typically  the  mathematical  model  is  a  set  of 
differential  equations  obtained  by  the  application  of  physical 
laws  that  govern  the  system.  The  response  of  the  system  to  a 
given  input  may  then  be  obtained  by  solving  the  differential 
equations. 

Development  of  the  equations  of  motion  begins  with  an 
application  of  Newton's  second  law  to  a  rigid  body  in  a  non¬ 
rotating  inertial-axis  system.  Motivation  is  given  to  make 
the  transformation  to  a  rotating,  body-axis  reference  system. 
Euler  angles  are  introduced  and  developed  in  some  detail  as  a 
means  of  referencing  body  axes  to  the  inertial  axes.  Finally 
a  brief  description  of  the  aerodynamic  model  is  presented  in 
the  form  of  non-dimensional  force  and  moment  equations. 

A.  APPLICATION  OF  PHYSICAL  LAWS 

Newton's  second  law  states  that  the  acceleration  of  an 
object  is  proportional  to  the  net  force  exerted  on  the  object. 
This  law  leads  directly  to  the  fundamental  concepts  of  the 
conservation  of  linear  momentum  and  the  conservation  of 
angular  momentum  [Ref.  20:p.  143].  In  a  non-rotating  inertial 
reference  frame  the  conservation  laws  can  be  expressed  in 
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terms  of  the  vector  equations: 

(1, 

(2) 

These  equations  state  that  the  summation  of  vector  forces  (F) 
acting  on  a  body  is  equal  to  the  time  rate  of  change  of  linear 
momentum  (mv)  and  the  summation  of  external  moments  (M)  acting 
on  a  body  is  equal  to  the  time  rate  of  change  of  angular 
momentum  (H) . 

The  externally-applied  aerodynamic  forces  and  moments 
acting  on  an  airframe  are  due  primarily  to  airflow  conditions 
and  control-surface  deflections.  The  force  and  moment 
equations  can  be  broken  down  into  vector  components  along  the 
longitudinal  (x)  ,  lateral  (y)  and  vertical  axis  (z)  of  the 
body  with  respective  velocity  components  u,v  and  w. 


d 

dt 

{mu) 

(3) 

L  = 

■4-^x 

dt  ^ 

II 

(4) 

Figure  1  illustrates  the  convention  for  the  forces  and  moments 
in  a  body-fixed  reference  frame  about  the  center  of  gravity. 
The  components  of  the  velocity  vector  V  are  u,  v  and  w  aligned 
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along  the  x,  y  and  z  axes  respectively.  The  moments  due  to 
aerodynamic  forces  are  L,  M  and  N. 


References  21  and  22  detail  the  development  of  the  moment 
equations  referred  to  a  moving  center  of  mass.  The  result  is 
a  vector  expression  for  angular  momentum  in  component  form 
given  by: 
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The  expression  for  the  angular  momentum  may  then  be  written 
as: 

(6) 

where  the  vector  n  contains  the  components  of  angular  velocity 
p,  q  and  r.  The  matrix  [I„]  is  a  symmetric  matrix  of  the 
moments  and  products  of  inertia  which  can  be  expressed  by: 
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with  components: 


j  iy^  +  z^)  dm 
(x^  +  z^)  dm 
I,=  f  (x^  +y2)  dm 


J  V 

ly^=fjyz)  dm 


(8) 


with  the  inertia  matrix  and  angular  velocity  vector 
defined,  equation  (2)  can  be  rewritten  as: 

Equation  (9)  indicates  that  in  a  fixed  reference  frame,  if  the 
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airframe  rotates  then  the  matrix  will  vary  with  time.  By 
making  the  transformation  to  a  rotating  coordinate  system  that 
is  fixed  to  the  airframe,  the  inertia  matrix  becomes  time 
invariant  and  the  analysis  becomes  much  simpler  [Ref  21]. 

In  making  the  transformation  from  the  inertial  reference 
to  the  rotating  body  reference,  the  vector  equations  for  P  and 
M  can  be  expressed  as: 

d 

F  =  m — (v)+ax(mv)  (10) 

dt 

M=—{H)^axH  (11) 

dt 


The  above  expressions  in  terms  of  the  body-axis  system  assume 
a  constant  mass  since  the  mass  of  the  aircraft  is  not  expected 
to  change  significantly  in  the  short  time  during  which  data 
are  collected.  The  components  of  the  vector  equations  (10) 
and  (11)  are  written  as; 

F^=m  (li  +qw-rv) 

Fy  =  m  (v+ru -pw)  (12) 

F^  =m  (w+pv-gu) 
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L=pj^-gj^y-rl„+gr 

M  =  -pj,/gjy-rjy, -^rp  ( J^-J^ )  +  (p2  -r2)  (13) 


B.  EULER  ANGLES 

The  development  in  the  previous  section  led  to  a  set  of 
six  equations  of  motion  derived  for  a  reference  system  that  is 
fixed  to  the  airframe.  Since  the  position  and  orientation  of 
the  airframe  cannot  be  described  relative  to  its  moving 
reference  frame,  it  is  necessary  to  develop  a  way  of  defining 
the  aircraft  attitude  relative  to  the  earth.  The  angles  that 
define  the  rotations  that  transform  the  orientation  of  the 
body  frame  to  the  earth-fixed  inertial  frame  are  called  the 
Euler  angles  [Ref.  22 :p.  88]. 

The  transformation  begins  by  aligning  the  body-axis  system 
with  the  inertial  system  (Figure  2)  .  Looking  down  the  Z  axis 
into  the  X,Y^  plane,  a  vector  R  in  that  plane  can  be  resolved 
into  X,  and  y,  components. 

R=x^jt  ^  y,9  (14) 

Rotate  an  angle  7  about  the  Z  axis  to  define  a  new  coordinate 
plane  X2Y2*  same  resultant  vector  R  can  be  expressed  in 
terms  of  the  components  of  the  new  coordinate  system. 
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J?=X2*  + 


(15) 


where 


Xj  =x,cosY+y^sinT 

y2=y,cos7-x,sinT  (16) 


In  matrix  form  the  new  coordinate  system  can  be  expressed  as 
a  product  of  the  transformation  matrix  T*  and  the  vector 
containing  components  x^,  Y\ 


^2 

cosT  sinT  0 

= 

-sinT  cosY  0 

yi 

^2 

0  0  1 

.^1. 

With  the  coordinate  system  now  aligned  with  the 
coordinate  system,  the  same  procedure  is  applied  by  looking 
down  the  Y2  axis  into  the  X2Z2  plane  and  observing  the  R  vector 
in  terms  of  components  X2  and  Zj  .  Rotate  an  angle  6  about  the 
Yj  axis  to  form  a  new  coordinate  system  Xj,  Yj,  Zj.  The  vector 
R  can  now  be  expressed  in  terms  of  components  X3,  y^,  Z3  by: 

Jl=X3±  +  Z3^  (18) 
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Figure  2  Euler  Angle  Rotation  About  the  Z  axis 

where  in  matrix  form  the  new  vector  is  a  product  of  the 
transformation  matrix  Tg  and  the  vector  in  the  Xj,  yj,  Zj 
system. 

^3  [co30  0  -sin0]  ^2 

y,  «  0  1  0  yj  (19) 

Zj  [sin0  0  COS0  Zj 

A  similar  procedure  is  applied  to  the  Xj,  Z3  coordinate 
system  by  rotating  an  angle  9  about  the  X3  axis.  The 


transformation  matrix  T,  is  found  to  be: 


10  0 


0  cos$  sin$ 
0  -sin®  cos® 


(20) 


A  complete  transformation  from  the  body  axes  to  the 
inertial  axes  by  three  successive  rotations  about  the  Euler 
angles  7,  8  and  #  leads  to  the  expression: 

=  (21) 

where  is  a  vector  of  components  in  the  body-axes  system  and 
X,-  is  a  vector  of  components  in  the  inertial-axes  system. 
Figure  3  shows  the  earth-fixed  axis  relative  to  the  body-  axis 
system  with  the  three  successive  rotations  about  the  angles  7, 
e  and  f. 

Using  the  relationship  in  equation  3.21,  the  inertial 
velocities  can  be  expressed  in  terms  of  the  Euler  angles  and 
the  components  of  velocity  in  the  body-axis  frame.  The 
shorthand  notation  of  References  21  and  22  where  represents 
Cos  6  gives: 


(22) 
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Figure  3  Euler  Angle  Orientation 

The  angular  rotation  rate  vector  n  has  body-axis 
components  p,  q  and  r  where  p  is  the  roll  rate  about 
the  longitudinal  axis,  q  is  the  pitch  rate  about  the  lateral 
axis  and  r  is  the  yaw  rate  about  the  vertical  axis.  Using  the 
relationships  of  the  transformation  matrices  T,,  and  T, 
given  above,  the  body-axis  components  of  the  angular  rotation 
rate  vector  n  can  be  expressed  in  terms  of  the  Euler-angle 
rates  by: 


p 

1  0  -Se' 

Q 

= 

0  C®  ^*^6 

© 

r 

0  ~S^ 

(23) 


Inverting  equation  (23)  gives  the  Euler-*angle  rates  in  terms 
of  the  body-axis  rotation  rates: 


1  Sin^bTanB  Cos^TanSi 

P 

© 

= 

0  CosO  -Sin® 

Q 

0  Sin^SecQ  Cos®  Sec© 

r 

(24) 


C.  FORCE  AND  MOMENT  EQUATIONS 

The  applied  forces  on  the  aircraft  can  be  expressed  in 
terms  of  aerodynamic,  gravitational  and  thrust  components. 
The  forces  acting  on  the  body  of  the  aircraft  relative  to  the 
earth  are: 

Fjf  =  -  mgSilid  +  T 

Fy  =  qsCy*mgSin<bCo^  (25) 

Fjj  =  gsCg+ mgCos<b  Cos© 

where  q-bar  is  the  dynamic  pressure,  s  is  the  reference  area, 
C^,  Cy,  Cj  are  the  force  coefficients  relative  to  the  aircraft 
state  and  T  is  thrust. 

The  gyroscopic  moment  due  to  rotating  parts  in  the  engine 
is  ignored  and  the  aerodynamic  moments  are  given  by: 
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M^  =  qscC^ 

My^gSCC„  (26) 

M^  =  gscC„ 

where  C^,  and  are  non-dimensional  coefficients  of  roll, 
pitch  and  yaw  moments. 

D.  STABILITY  AND  WIND  AXES 

For  stability  and  control  analysis  it  is  more  convenient 
to  consider  the  velocity  of  the  aircraft  relative  to  the  air 
[Ref  18 :p  8].  Also  it  is  more  practical  to  express  equations 
in  terms  of  angle  of  attack,  a  and  sideslip  angle,  3  since 
these  angles  are  measured  directly  during  wind  tunnel  testing 
and  flight  testing.  Angle  of  attack  and  sideslip  are  defined 
by; 

o*tan-i^  P=sin-^-^  (27) 

where  V  is  the  total  wind-relative  velocity  with  magnitude 
given  by; 

V=}/u^  +  v^  +  w^  (28) 


From  Figure  4  the  component  velocities  can  be  expressed  in 
terms  of  the  angle  rotations  through  alpha  and  beta  and  the 
total  wind-relative  velocity  V  by  the  following  process. 
Rotate  about  the  w  axis  by  an  angle  3; 
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(29) 


^2 

Cp  -5P 

0 

^2 

= 

sp  cp 

0 

’^2. 

0  0 

1 

Rotate  down  about  the  v  axis  by  an  angle  a: 


“3 

Ca  0  -Sa 

“2 

^3 

= 

0  10 

Sa  0  Ca 

*^2 

(30) 


The  component  velocities 


can  then  be  expressed  by: 


u 

V 

V 

hT 

II 

0 

w 

0 

(31) 


Figure  4  Flow  Angles  Relative  to  Alpha  and  Beta 
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Through  a  similar  application  of  the  angle  relationships 
the  axial  forces  can  be  expressed  in  terms  of  the  non- 
dimensional  coefficients  Cp  and  and  related  to  the  stability 
axes  and  the  wind  axes  by  a  transformation  matrix. 


So. 

~Ca 

-Ca 

-Sa, 

c 

^wtnd 


cp  -sp 

5P  Cp 

Cy 

(33) 


E.  BODY  AXIS  EQUATIONS  OF  MOTION 

Constructing  the  body-axis,  six-degree-of-freedom 
equations  of  motion  begins  by  differentiating  equations  (27) 
and  (28)  to  obtain  expressions  for  the  time  derivatives  of  V, 
alpha  and  beta.  Equating  equations  (12)  and  (25)  yields: 

u=  -2®  r  -qrv  +  rv-gsin0  +  — 
m  ^  m 

v=  —  C  -ru  +  pw+gsin®cos0  (34) 

w=  Cj-pv  +  gu  +  gcos4>cos0 


Substitute  equation  (31)  for  u,  v,  and  w  and  equation  (34)  for 
the  time  derivatives  of  u,  v,  and  w.  Use  equations  (32)  and 
(33)  to  substitute  in  for  the  force  coefficients  Fx,  Fy  and 
Fz.  The  construction  is  completed  by  substituting  equation 
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(33)  to  substitute  in  for  the  force  coefficients  Fx,  Fy  and 
Fz.  The  construction  is  completed  by  substituting  equation 
(26)  into  equation  (13).  The  resulting  nine  equations  can 
then  be  expressed  as: 


V=  +  g  (cosOcosGsinacosP  -t-sin^cosBsinP 

T 

-sin0cosacosP )  -  — cosacosp 

in 


(35) 


d  =  -- 


gs 


(tanP)  (pcosa  +rsina) 


/nVcosp 

+  — -2-— (cos®cos0cosa +  sin6sina) 

VtosP  mV^cosP 


(36) 


P  =  -2^Cy  +psina -rcosa +  ^cosPsin®cos0 

(gcosasin0 'gsinacos<>cos0  +  — coso) 
V  m 


(37) 


<^I^-tJ^^=qsbC^*qr(Iy-I^)  +{q^-x^)  ly^ 

(38) 

^pql^-rpl^ 

(39) 

^qrl^-pqly^ 

,-qIy^^rI^=qsbC„+pq{I^-Iy)  +  (p^-q^) 

(40) 

*rply^-qzl„ 

^  =p  +  gtan0sin®  +  rtan0cos® 

(41) 

0  =  gcos®  -rsinC> 

(42) 
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=rcos®sec0  +  gsin®sec0 


(43) 


The  nine  equations  listed  above  are  valid  for  a  rigid 
vehicle  in  a  constant  wind  with  the  following  assumptions. 

1)  Flat,  non-rotating  earth 

2 )  Constant  mass 

3)  Thrust  vector  along  the  X  axis 

4)  Body  axis  referenced  to  the  aircraft  center  of  gravity 


Singularities  occur  at  e  =  ±  90  degrees,  6  =  ±  90  degrees  and 
at  zero  velocity. 


F.  AERODYMAMIC  MODEL 

Equations  (12)  and  (25)  can  be  combined  to  give  a 
representation  of  the  forces  acting  on  the  aircraft.  Changes 
in  the  forces  from  a  given  reference  condition  can  be 
represented  by  a  Taylor  series  expansion  as  a  function  of  the 
motion  variables.  In  the  development  of  the  aerodynamic 
model,  it  is  assumed  that  any  change  in  the  longitudinal  force 
is  a  function  of  angle  of  attack,  pitch  rate  and  longitudinal 
control.  Changes  in  the  longitudinal  forces  and  moment  can 
then  be  represented  by; 


du  aw  00^ 


A6. 


du  dw  dw 

AM= A  u  + -3^  A  w+ —  A 
du  dw  dw 


dM  .  du  y.  • 


(44) 
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Likewise  it  is  assumed  that  the  lateral-directional 
aerodynamic  forces  and  moments  are  functions  only  of  sideslip 
angle,  roll  rate,  yaw  rate  and  lateral-directional  control- 
surface  changes.  Expansions  for  AY,  AL  and  AN  can  then  be 
expressed  in  terms  of  the  motion  variables  Av,  Ap,  Ap,  Ar  and 
A6  by; 


dv  dp  dr  do^  ^ 

at  dL  dL  dL  Si-  A  X  ^  SL  A  X 

8^  ^ *  a?  * -55;  "  is; 

AAr=-^Av+-^Ap  +  -^Ar  +  -^  A6  +-^AA 
dv  dp  ^  dr  d6^  ^  36,  ' 


(45) 


In  equations  (44)  and  (45)  a  conventional  three-axis 
control  system  consisting  of  aileron,  rudder  and  elevator  is 
assumed.  Additional  control  terms  may  be  added  for  non- 
conventional  configurations.  Reference  23,  Chapter  IV 
clarifies  the  significance  of  each  of  the  derivative  terms  in 
the  expansions  and  explains  the  physical  significance  by 
relating  the  longitudinal  and  lateral-directional  forces  and 
moments  to  the  aerodynamic  derivatives  and  motion  variables. 

Stability  and  control  derivatives  in  dimensional  form  are 
convenient  for  application  of  modern  control  theory  methods 
when  solving  flight  dynamic  problems.  References  21  and  22 
use  the  dimensional  form  in  the  development  of  the  state-space 
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model.  The  non-dimensional  coefficients,  however,  provide 
a  convenient  means  of  accounting  for  aircraft  speed  and  size, 
and  for  air  density  since  they  are  insensitive  to  dynamic 
pressure  over  large  ranges  [Ref.  24:p.  113].  Table  1 
summarizes  the  longitudinal  stability  and  control  derivatives 
in  terms  of  non-dimensional  coefficients  and  Table  2 
summarizes  the  lateral-directional  stability  and  control 
derivatives.  The  derivatives  represent  the  change  in  an 
aerodynamic  force  or  moment  coefficient  due  to  a  change  in 
state  or  control  input. 
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TABLE  1  LONGITUDINAL  STATES  AND  DERIVATIVES 


Longitudinal  Stability  and  Control  Derivatives 

State  or  Control 

Normal  Force 

Pitch  Moment 

Angle  of  Attack 

CL. 

Cm. 

Pitch  Rate 

CLq 

Cm^ 

Elevator  Control 

CL,. 

C^ie 

TABLE  2  LATERAL-DIRECTIONAL  STATES  AND  DERIVATIVES 


Lateral-Directional  Derivatives 


State  or  Control 

Side  Force 

Roll  Moment 

Yaw  Moment 

Sideslip 

CYe 

CIb 

CUq 

Roll  Rate 

CYp 

Cnp 

Yaw  Rate 

CY, 

Clr 

CHp 

Aileron  Control 

CY„ 

cij. 

Cn,. 

Rudder  Control 

CY,, 

Cl,r 

Cn,p 
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IV.  MAXIMUM  LIKELIHOOD  ESTIMATION 


A.  ESTIMATION  THEORY 

The  problem  of  obtaining  unknovm  aircraft  stability  and 
control  derivatives  through  parameter  estimation  techniques  is 
theoretically  straight  forward.  In  the  previous  chapter  it 
was  emphasized  that  the  most  important  part  of  this  analysis 
is  deriving  a  reasonable  mathematical  model  that  consists  of 
a  set  of  dynamic  equations  that  contains  the  unknown 
parameters.  The  system  can  then  be  excited  by  a  pre> 
determined  input  specifically  designed  so  that  the  response  of 
the  system  will  allow  the  values  of  the  parameters  to  be 
determined.  The  system  response  to  the  pre-determined  input 
is  measured.  Values  of  the  parameters  are  then  determined  by 
the  constraint  that  the  computed  response  of  the  model  match 
the  measured  response  of  the  system.  This  deterministic 
approach  to  parameter  identification  is  complicated  by  the 
presence  of  other  elements  such  as  process  noise  and 
measurement  noise  that  must  be  considered  in  real-world 
systems . 

The  mathematical  model  derived  in  the  previous  chapter  is 
based  on  the  physical  laws  that  govern  the  system  and  is 
considered  to  be  well  defined.  In  deriving  any  model  of  a 
physical  system  it  is  necessary  to  make  some  compromise 
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between  simplicity  and  accuracy.  A  model  that  is  too 

complicated  will  require  a  higher  expenditure  of  resources  and 

a  model  that  is  too  simple  will  not  have  sufficient  accuracy 

to  be  useful.  The  objective  is  to  define  a  mathematical  model 

that  represents  the  essential  aspects  of  the  physical  system 

and  will  provide  reasonably  accurate  predictions  of  the 

dynamic  behavior  [Ref.  21:p.  67].  The  term  "reasonably 

accurate"  Implies  that  modeling  errors  are  inevitable.  In 

Reference  17  Maine  and  Iliff  state: 

There  is  no  comprehensive  theory  of  modeling  error 
available.  Any  modeling  error  is  simply  treated  as  state 
and/or  measurement  noise,  in  spite  of  the  fact  that  the 
modeling  error  may  be  deterministic  rather  than  random. 
This  procedure  has  not  been  rigorously  justified,  but 
combined  with  careful  choice  of  the  model,  is  probably  the 
best  approach  available. 

These  unavoidable  modeling  errors  are  intentionally  accepted 
in  an  effort  to  produce  a  simplified  but  practically  useful 
simulation. 

The  presence  of  both  measurement  noise  and  state  noise 
changes  the  nature  of  the  problem  from  a  deterministic  process 
to  a  stochastic  process.  Measurement  noise  arises  from  the 
realization  that  the  system  response  cannot  be  perfectly 
measured  and  sensor  inaccuracies  are  certain.  State  noise  is 
random  excitation  of  the  system  from  external,  unmeasured 
sources.  The  noise  is  a  disturbance  that  acts  uninterruptedly 
upon  the  system  making  it  impossible  to  predict  exactly  how 
the  system  will  behave.  Consequently  the  unknown  parameters 
must  be  estimated  by  some  statistical  criterion  [Ref  17 :p  8]. 
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The  real  value  of  the  MMLE3  program  is  that  the  algorithm 
accounts  for  both  state  and  measurement  noise. 


B.  MODIFIED  MAXIMUM  LIKELIHOOD  VERSION  THREE 
1.  General  System  Model 

In  the  Modified  Maximum  Likelihood  Estimation  (MMLES) 
method  it  is  assumed  that  for  the  correct  values  of  the 
unknown  parameters,  the  system  is  correctly  described  by  the 
dynamic  model  developed  in  the  previous  chapter.  The  model 
can  be  expressed  in  state-space  format  as  a  function  of  the 
states,  inputs  and  unknown  parameters  represented  by  z,  u  and 
^  respectively.  The  discretized  model  in  state-space  form  is: 

+  (46) 

where  the  plant  matrix  A  and  the  control  matrix  B  contain 
values  of  the  unknown  parameters  and  u  is  the  state  noise 
vector.  The  initial  condition  is  given  by: 


Measurements  are  assumed  to  be  made  at  discrete  time 
intervals.  The  measurement  equation  is: 


(48) 
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The  matrix 

C  is 

the 

measurement  matrix  and  v 

k  is 

the 

measurement 

noise 

vector . 

The  state 

noise,  e. 

and 

the 

measurement 

noise, 

are 

assumed  to 

be  zero 

mean 

and 

uncorrelated  [Ref.  25:p.  3].  In  the  presence  of  white  noise 
the  estimation  is  a  stochastic  process.  For  each  estimate  of 
the  unknown  parameters  in  the  vector  there  is  associated 
with  the  estimate  a  probability  that  the  values  in  the 
response  history  are  close  to  the  measured  values.  The 
maximum  likelihood  estimates  are  the  estimates  that  maximize 
that  probability. 

2 .  Kalman  Estimator 

A  Kalman  estimator  is  used  to  provide  estimates  for 
the  states  while  minimizing  the  effects  of  the  process  and 
measurement  noise.  The  difference  between  the  measured  value 
and  the  Kalman  estimate  at  a  given  time,  t(i),  is  called  the 
residual  and  is  given  by: 


(49) 


The  subscript  i  indicates  that  the  estimate  is  a  function  of 
the  unknown  parameters.  To  maximize  the  probability  that  the 
computed  time  history  is  close  to  the  measured  time  history  it 
is  necessary  to  minimize  the  criterion  given  by  [Ref  25:p  3]: 


(50) 
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Where  R  is  the  innovations  co-variance  matrix.  From  equation 
(50)  it  can  be  seen  that  the  effect  of  the  state  and 
measurement  noise  is  only  to  increase  the  value  associated 
with  the  criterion  J({). 

The  discrete,  stochastic  estimation  is  a  recursive 
process  that  makes  an  estimate  after  each  measurement.  The 
new  estimate  depends  upon  the  prior  estimate  and  the  new 
measurement.  In  terms  of  a  Kalman  estimator  equation,  the 
prediction  and  correction  steps  can  be  written  as: 


^k=^k*K[Z^-Z^] 


(51) 


Equation  (51)  states  simply  that  the  corrected  value  is  equal 
to  the  old  estimate  plus  a  confidence  factor  times  the 
residual.  The  gain  or  confidence  factor  is  chosen  by  solution 
of  the  Riccati  equation  such  that  the  errors  associated  with 
the  states  are  driven  towards  the  convergence  criterion  as 
quickly  as  possible.  MMLE3  uses  only  the  steady  state 
solution  of  the  Riccati  equation.  A  complete  explanation  of 
the  discrete,  stochastic  process  can  be  found  in  Chapter  6  of 
Reference  27. 

3 .  Cramer-Rao  Bounds 

In  any  stochastic  process  some  measure  of  accuracy 
is  necessary.  In  the  MMLE3  estimator  the  Cramer-Rao  bound  is 
a  measure  of  the  reliability  of  the  estimate  based  on  the 
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information  obtained  from  the  dynamic  maneuver.  In  simpler 
terms,  the  Cramer-Rao  bound  is  an  estimate  of  the  standard 
deviation  of  the  estimates.  The  measure  of  reliability  is 
improved  when  a  sufficient  number  of  measurements  are  provided 
and  the  assumptions  used  in  the  prediction  model  are  valid. 
The  Cramer-Rao  bounds  are  the  square  roots  of  the  diagonal 
elements  of  the  information  matrix.  A  precise  mathematical 
description  of  the  information  matrix  and  the  calculation  of 
the  Cramer-Rao  bound  can  be  found  in  References  25  and  27. 
For  the  purposes  of  this  study  it  is  enough  to  know  that  the 
Cramer-Rao  bound  is  a  good  indicator  of  the  accuracy  of  an 
estimated  parameter. 

4 .  A-Priori  Information 

Wind  tunnel  data  from  the  study  conducted  by  the 
NASA  Langley  Research  Center,  Flight  Dynamics  Branch,  provide 
some  insight  into  the  aerodynamic  characteristics  of  the  BQM- 
14 7A.  Measurements  tabulated  in  Reference  2  are  useful  in 
providing  estimates  for  some  of  the  stability  and  control 
parameters.  The  MMLE3  parameter  estimation  routine  allows  the 
estimation  algorithm  to  consider  this  a-priori  information  by 
adding  a  quadratic  penalty  function  to  the  criterion  [Ref. 
27 :p.  SR-4].  Equation  (50)  can  then  be  re-expressed  by; 
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(52) 


N 

2  i.i 

i(5-5o)' W  a-io) 

2 

where  5q  is  the  vector  of  a-priori  values  and  W  is  the 
weighting  matrix.  The  a-priori  information  can  be  used  to 
speed  the  convergence;  however,  the  results  will  be  biased 
towards  the  weighted  estimates  [Ref  27 :p  ST-8].  A  good 
technique  to  avoid  the  problems  associated  with  biased 
estimates  is  to  remove  the  weighting  after  initial 
convergence.  The  initial  weighted  values  should  be  close 
enough  to  the  maximum  likelihood  estimates  to  allow  rapid 
convergence  to  an  unbiased  estimate  when  the  weighting  is 
removed  and  additional  iterations  are  allowed. 

5.  Linearized  Equations  of  Motion 

In  the  analysis  of  aircraft  stability  and  control 
it  is  often  more  practical  to  restrict  the  area  of  interest  to 
steady-state  motion  and  the  response  of  the  aircraft  to  small 
perturbations  about  a  steady-state  condition  [Ref  24  :p  33]. 
Equations  (35)  through  (43)  in  Chapter  3  are  a  full  set  of 
coupled,  non-linear  equations  that  completely  describe  the 
dynamic  behavior  of  the  aircraft.  The  simplification  of  the 
model  by  applying  small-angle  theory  significantly  reduces  the 
complexity  of  the  model  without  necessarily  compromising 
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accuracy.  For  a  wide  range  of  practical  applications  on  well- 
behaved  systems,  a  linear  model  is  sufficient. 

References  22  through  25  provide  details  in  applying 
small-angle  theory  to  the  non-linear  equations  of  motion  from 
which  the  linearized  equations  used  in  MMLE3  are  derived. 
From  equations  (36) ,  (39)  and  (42) ,  a  small  angle 
approximation  leads  to  the  simplified  longitudinal  equations: 


A*- 


^  V 


(53) 


“  «# 


(54) 


(55) 


In  the  approximation,  a  plane  of  symmetry  about  the  X-Z  plane 
is  assumed,  making  ly^  and  equal  to  zero.  The  aerodynamic 
coefficients  in  the  above  equations  are: 


=  -2£+r 

mm,  m^  "‘q  2V 


(56) 


(57) 


From  equations  (56)  and  (57)  the  vector  5  is  expressed  as  a 
function  of  the  unknown  stability  and  control  derivatives  by: 


5  ^L»* 


(58) 


From  the  application  of  the  small  angle  approximation  to 
equations  (37) ,  (38) ,  (40)  and  (41) ,  the  simplified  lateral- 
directional  equations  are: 

=  gsjbCj  -*•  gr  ( ly-I^)  ( 60 ) 

^^2  -P^xz  =  QsbC„  *pq  ( -  ly)  -grJ^^  ( 61 ) 


4>®p  +  rcos<J>tan0  +  gsin4>tan0 


(62) 


The  aerodynamic  coefficients  for  the  lateral-directional 
equations  are: 


=  C„  0  +  C.,  -*•  6  +  C 

Y  rtf'  ^i>2V  2V  ' 


(63) 


C,  =  C,  6  +  C,  +  C,  —  +  C,  6  +  C, 


(64) 
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(65) 


rb 


■■C„  B  +C„  ^  +C„  —  ■••C- 6  +C_ 
Up  2^  2V  ^  ^ 


The  "delta"  parameter  is  a  general  expression  for  control 
input  and  must  be  expanded  to  consider  all  inputs  that 
contribute  to  the  total  force  or  moment  coefficient.  The 
contributions  to  the  total  side  force  coefficient  from  roll 
rate  and  yaw  rate  are  small  compared  to  the  other  contributing 
parameters  and  are  usually  ignored.  The  vector  of  unknown 
parameters,  can  be  expressed  as  a  function  of  the  lateral- 
directional  stability  and  control  derivatives  by: 


5=5  (Cy.  C, 


"p 


I  •r 


(66) 


C.  pBst 

For  many  aircraft  applications  the  linearized 
equations  of  motion  provide  a  reasonably  accurate  model  of  the 
system  for  small  disturbances.  For  some  applications  however, 
the  linear  model  is  not  sufficient  and  a  more  accurate  non¬ 
linear  model  is  required.  Flight  testing  at  the  Dryden  Flight 
Research  Facility  highlighted  the  limitations  of  the  linear 
model  when  applied  to  maneuvers  at  extreme  flight  conditions. 
In  1987  Murray  and  Maine  [Ref  9]  introduced  the  pEst  program 
for  parameter  estimation  in  non-linear  dynamic  systems. 


3 


1.  General  System  Model 


Conceptually,  parameter  estimation  in  non-linear 
dynamic  systems  is  the  same  as  the  parameter  estimation  in 
linear  systems.  The  general  system  model  is  separated  into  a 
continuous-time  state  equation  and  a  discrete  response 
equation  represented  by  [Ref.  9:p.  3]: 

i(C)  =f[x{t)  ,u{t)  ,i]  (67) 

z  (c^)  =g[x(ti) ,  uit^)  ,  1]  (68) 

The  cost  function  is  the  criterion  by  which  the  model's 
computed  response  is  measured  against  the  measured  time 
history  response.  The  pEst  program  defines  quantitatively  the 
criterion  by: 

Jii)  E  W  [ziti)  -z(ti)]  (69) 


Where  n^  and  n^  are  the  number  of  response  variables  and  the 
number  of  measured  time  history  points. 

2 .  Estimator 

The  feedback  process,  like  that  used  in  the  MHLE3 
estimator,  is  recursive.  The  non-linear  equations  used  in 
pEst  preclude  using  the  discrete-time  Kalman  estimator 
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formulation  of  MMLE3;  however,  the  feedback  term  is  still 
proportional  to  the  difference  between  the  measured  response 
and  the  computed  response  [Ref.  9:p.  3]: 

=x(t^)  *k[z(t^)  - (70) 

3 .  Parameters 

The  parameters  are  the  values  that  can  be  estimated  by 
pEst.  Each  of  the  parameters  has  five  attributes:  current 
value,  estimation  status,  predicted  value,  Cramer-Rao  bound 
and  change  in  value  from  the  previous  iteration  [Ref  9:p  11]. 
The  initial  current  value  is  a  best  estimate  of  the  parameter 
entered  by  the  user.  Typically  the  initial  current  value 
input  is  derived  from  a  source  such  as  DATCOM,  [Ref.  30  ]  or 
it  may  be  the  result  of  measured  data  from  either  wind  tunnel 
studies  or  flight  tests.  The  estimation  status  indicates  to 
the  user  which  parameters  are  being  estimated.  A  parameter  is 
either  active  or  inactive.  The  predicted  value  is  the  value 
of  the  parameter  computed  by  pEst  and  it  is  used  to  update  the 
current  value  for  each  successive  iteration.  The  Cramer-Rao 
bound  and  the  change  in  parameter  value  are  defined  in  the 
estimation  process.  The  Cramer-Rao  bound  is  a  measure  of  the 
accuracy  of  the  estimate.  A  list  of  the  names  and  a 
description  of  the  parameters  available  in  pEst  can  be  found 
in  Section  B.16,  Reference  9.  For  this  study,  the  parameter 
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vector,  is  the  same  as  the  parameter  vector  used  in  the 
MMLE3  routine. 

4.  States  end  Responses 

The  attributes  o£  the  states  are  status  and 
integration  limit.  The  attributes  of  the  response  variables 
are  status  and  weighting.  If  a  state  is  active,  the  state 
equation  is  integrated  and  its  value  is  used  in  the  equations 
of  motion.  If  the  value  of  the  state  exceeds  the  integration 
limit,  the  integration  is  terminated.  If  the  response 
variable,  or  output  variable,  is  active  then  its  time  history 
values  are  computed  and  placed  in  a  special  output  file.  The 
weighting  of  a  response  variable  specifies  the  value  to  be 
used  in  the  W  matrix  of  the  cost  function.  Weighting  the 
response  variables  based  on  a-priori  knowledge  may  help  the 
routine  converge  more  rapidly;  however,  the  results  may  be 
biased  towards  the  weighted  estimates.  Table  4.13  of 
Reference  24  shows  the  relative  importance  and  prediction 
accuracies  of  the  stability  derivatives  using  theoretical 
methods.  The  table  may  be  used  as  a  guideline  if  weighting  is 
desired. 

5.  Equations  of  Motion 

Equations  (35)  through  (43)  define  the  non-linear 
equations  of  motion  used  in  the  pEst  routine.  The  assumption 
of  symmetry  about  the  X-Z  plane  that  was  used  in  MMLE3  is  not 
used  in  pEst;  however,  the  assumptions  of  fixed  aircraft 
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geometry  and  constant  mass  are  valid  [Ref.  9:p.  17].  In  this 
study,  equations  (36)  through  (42)  are  integrated  to  obtain 
the  predicted  state  variables  in  the  vector  x-tilde.  A  list 
of  the  14  response  equations  can  be  found  in  Section  6.1.4  of 
Reference  9. 
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V.  MMLE3  APPLICATION 

Recent  thesis  work  by  LCDR  Robert  Graham  [Ref.  6] 
incorporated  a  PC-based  parameter  estimation  capability  into 
the  NPS  research  program.  The  Maine  and  Iliff  MMLE3  code  has 
been  used  for  a  number  of  years  for  large  scale  work  where  the 
use  of  mainframe  computer  systems  was  justified  [Ref.  9:p.  1] . 
Research  by  Erickson  [Ref.  8]  showed  that  for  less  demanding 
work,  a  PC-based  parameter  estimation  routine  for  a 
longitudinal  model  generally  provides  comparable  results  when 
compared  to  more  complex  mainframe  programs.  The  code  written 
by  Graham  uses  the  commercially  available  software  package, 
MATLAB*,  to  build  interactive  files  for  both  the  longitudinal 
and  lateral-directional  models.  The  MATLAB  State-Space 
Identification  Tool,  (SSID)  [Ref.  28],  comprises  a  main 
Modified  Maximum  Likelihood  Estimator  (MMLE)  program,  along 
with  several  supporting  functions  for  parameter  estimation. 
The  supporting  functions  are  used  to  convert  the  parameter 
vector  into  the  state-space  model  and  to  arrange  the  necessary 
inputs  for  the  MMLE  program.  A  complete  listing  of  the  code 
written  by  Graham  can  be  found  in  Appendices  E  and  F  of 
Reference  6.  The  modifications  to  the  supporting  functions  in 
Reference  6  required  for  this  application  are  discussed  in  the 
following  sections.  Figure  5  [Ref.  9:p.  43],  shows  how  the 
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various  MATLAB  files  interact  in  the  parameter  estimation 
routine. 


Figure  5  Interactive  MATLAB  .M  Files 


A.  DATA  ACQUISITION 
1.  Flight  Test 

The  data  made  available  for  this  study  were  provided 
by  NASA  Langley  Research  Center  and  were  the  result  of  flight 
tests  conducted  at  the  NASA  Langley  Plumtree  Test  Site  in 
March,  1992.  This  flight  test  comprised  a  single  set  of 
flight  maneuvers  consisting  of  pitch,  roll  and  yaw  doublets. 

A  partial  sequence  of  events  is  listed  in  Table  3.  The  roll 
and  yaw  doublets  were  merged  into  a  single  input  file  so  that 
all  of  the  lateral-directional  modes  would  be  excited. 
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Figures  6  and  7  show  a  typical  control  input  for  the  doublet 
maneuvers . 

2 .  Atmospheric  Data 

Specific  atmospheric  data  for  the  date  and  time  of  the 
flight  tests  were  unknown.  Best  estimates  of  atmospheric 
conditions  for  the  test  site  were  from  a  local  area  weather 
briefing  and  were  recorded  on  the  flight  card  by  the  flight 
test  pilot. 

3 .  Sensor  Placement 

Twin  alpha  and  beta  sensors  were  symmetrically  located 
about  the  longitudinal  axis  and  2.5  feet  in  front  of  the 
vehicle  center  of  gravity  (eg) .  An  averaging  scheme  was  used 
on  the  two  sets  of  measured  data  and  a  single  set  of  data  for 
alpha  and  beta  was  provided.  Accelerometers  for  X,  Y  and  Z 
axes  were  located  on  the  eg  point.  A  rate  sensing  gyro  was 
located  on  the  waterline,  aft  of  the  eg.  Details  on  sensor 
characteristics  were  not  available. 

The  aircraft  equations  of  motion  developed  in  Chapter  III 
assume  that  the  sensors  are  located  at  the  eg.  Since  several 
of  the  observed  signals  are  functions  of  the  position  of  the 
sensors,  it  is  necessary  to  correct  the  sensor  measurement  for 
the  offset.  For  the  longitudinal  model  the  measured  variables 
and  A„  are  corrected  to  the  eg  by  [Ref.  18 ;p.  26]; 


a  =a  + _ —Q 

ew  «  ..  ^ 


(71) 
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(72) 


=  A,- 


g  9 


For  this  application,  is  positive  forward  of  the 
eg  and  was  given  a  fixed  value  of  +2.5  feet.  Since  the  X-and- 
Z  axis  accelerometers  were  located  at  the  eg,  the  value  of  X^ 
and  is  zero.  For  the  lateral-directional  model,  the 
measured  variables  13^  and  are  corrected  to  the  eg  by  [Ref 
18;p  27] ; 


(73) 


(74) 


The  X  distance  to  the  beta  probe  was  set  at  +2.5  feet  and  the 
value  of  X-„  and  Z-„  was  set  to  zero. 

4 .  Measured  signals 

A  total  of  twenty-four  signals  were  measured  and 
recorded  during  the  flight  test.  These  signals  are  presented 
in  Table  4.  Details  of  how  the  data  were  filtered,  digitized, 
time  tagged  and  recorded  were  not  available. 
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B.  USER  FUNCTION 

The  PC-based  MMLE3  application  requires  a  user-written 
function  to  convert  the  parameter  vector  into  a  state-space 
system  model  [Ref.  28 :p.  SR-3].  The  following  two  sub¬ 
sections  describe  the  creation  of  the  P2SS.m  file  for  this 
study. 

1.  Longitudinal  Modal 

In  the  general  system  model  given  by  equation  (46), 
the  response  of  the  system  is  a  function  of  the  state  vector, 
X,  and  the  control  input,  u.  Since  the  measured  flight  test 
data  do  not  include  pitch  angle  and  in  the  short-period 
approximation  forward  velocity  is  assximed  to  be  constant,  the 
simplified  longitudinal  model  can  be  reduce  to  equations  (53) 
and  (54)  for  a  and  q.  These  two  equations  contain  terms  which 
are  neither  functions  of  the  state  variables  or  the  control 
inputs.  Consequently,  the  control  input  was  augmented  with  a 
unity  input  to  incorporated  these  terms  in  the  state-space 
model.  The  state  vector  is  represented  by: 


The  control  input  is  represented  by: 


(76) 
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In  the  state-space  formulation  the  vector  p  is  the  vector 
of  stability  and  control  parameters.  This  vector  is  expressed 
in  terms  of  its  elements  by: 

(77) 


The  plant  matrix  can  now  be  expressed  in  terms  of  the  elements 
of  the  parameter  vector  by: 


- — P(l) 

gwVcosaj 

^f£p{2) 


(78) 


The  control  matrix  is  expressed  in  terms  of  the  elements  of 
the  parameter  vector  by: 


- 2® - p(4) 

gwVcosa^ 


-^P(5) 


QS 

gwVcosa^ 
(P(2)  Sj 


(p(l)  +p(4)  6^^)  -Qi 
■^P(5)  «e,)  <7i 


(79) 


In  the  longitudinal  application  it  is  desired  to  evaluate 
the  dependence  of  any  change  in  normal  acceleration  upon 
longitudinal  control  input.  The  equation  for  normal 
acceleration  in  terms  of  non-dimensional  stability  and  control 
parameters  is  derived  in  Chapter  4  of  Reference  22.  The 
general  system  output  equation  can  be  expressed  in  terms  of  a 
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parameters  is  derived  in  Chapter  4  of  Reference  22.  The 
general  system  output  equation  can  be  expressed  in  terms  of  a 
parameter  vector,  y,  that  contains  the  variable  A„  in  addition 
to  the  states  a  and  q.  The  output  and  measurement  equations 
of  the  form  given  in  equation  (48)  are  augmented  by  a  feed¬ 
through  matrix,  D.  The  output  matrix,  C,  can  be  expressed  by: 


C= 


1 

0 


QS 

mg 


p(l) 


o' 

1 

0 


(80) 


The  feed-through  matrix,  D,  is  given  by: 


0 

0 


-2^p(4) 

gw 


(81) 


The  subscripted  states,  and  q^,  are  the  elements  of  the 
initial  state  vector  Xq.  The  input  for  the  initial  vector 
comes  directly  from  the  first  element  in  the  measured  data 
file  for  each  of  the  states. 

The  longitudinal  model  given  above  specifies  the  state 
equation  in  continuous  form.  In  computing  the  continuous-to- 
discrete-time  conversion,  the  sample  time,  dt,  is  required. 
Since  MATLAB  will  not  allow  a  global  variable  to  be  used  in  a 
function  argument,  the  P2SS.m  file  must  be  modified  to  allow 
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2 •  Lataral-Directlonal  Modal 


In  the  lateral-directional  model  the  state  vector,  x, 
is  represented  by: 


jf=  [B  p  r  0]^ 


(82) 


The  control  input  vector,  u,  is  represented  by: 


u  = 


(83) 


The  subscripted  variables  and  5^  represent  the  aileron  and 
rudder  control  inputs.  The  input  vector  is  augmented  by  a 
unity  input  to  incorporate  the  additional  terms  that  influence 
lateral  acceleration  in  the  measurement  equation.  The 
measurement  vector  is: 

y=[Bpr0Ay,]^  (84) 


The  lateral-directional  parameter  vector,  p,  is: 


c 

a 


(85) 


The  only  change  to  the  NPSP2SS4.M  file  in  Reference  6 
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was  to  input  the  time  step,  dt,  for  the  sampling  rate  of  the 
measured  data  file. 

From  Table  4  it  can  be  seen  that  roll-angle  was  not 
measured  during  flight  test.  Since  the  spiral  mode  is 
dominated  by  the  roll  angle  component,  it  was  necessary  to 
numerically  integrate  the  roll  rate,  p,  in  order  to  provide  a 
roll  angle  input  estimate.  For  N  intervals  with  equal 
spacing,  the  extended  trapezoidal  rule  can  be  written  as: 


N-1 


f{a)  +252  f  (a+j  (dt) )  +f  (b) 
j-i 


+  E 


(86) 


A  complete  development  of  the  extended  trapezoidal  rule  and  a 
discussion  of  the  error  term,  E,  can  be  found  in  Chapter  IV  of 
Reference  29. 

C .  INITIALI ZATION 

The  initialization  macro  is  a  supporting  function  that 
imports  the  measured  time  history  file,  prompts  the  user  for 
the  physical  dimensions  of  the  aircraft  and  prompts  for  the 
initial  values  of  the  parameters  that  are  to  be  identified. 

1.  Physical  Characteristics 

Table  5  presents  the  aircraft  physical  characteristics 
that  were  input  for  this  study.  To  be  consistent  with  the 
physical  characteristics  used  in  the  NASA  static  wind-tunnel 
tests,  the  theoretical  root  chord  was  used  as  the  reference 
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1.  Physical  Characteristics 

Table  5  presents  the  aircraft  physical  characteristics 
that  were  input  for  this  study.  To  be  consistent  with  the 
physical  characteristics  used  in  the  NASA  static  wind-tunnel 
tests,  the  theoretical  root  chord  was  used  as  the  reference 
wing  chord.  Moments  of  inertia  and  sensor  offset  information 
were  provided  by  the  Marine  Corps  UAV  Project  Office, 
Intelligence  Systems,  Quantico,  VA. 

2.  Initial  Estimates 

The  SSID  Modified  Maximum  Likelihood  Estimator 
requires  a  set  of  predicted  aerodynamic  characteristics  to  be 
used  as  starting  values  for  the  iterative  algorithms.  In  this 
study  a  combination  of  wind  tunnel  data  and  analytic  estimates 
was  used  as  a  source  of  predicted  values.  Tables  6  and  7 
present  the  parameter,  the  predicted  value  of  the  parameter 
and  the  source  of  the  predicted  value  for  the  longitudinal  and 
lateral-directional  case  studies. 

D.  MODIFIED  MAXIMUM  LIKELIHOOD  ESTIMATOR 

The  Modified  Maximum  Likelihood  Estimator  is  the  main 

program  used  by  the  MATLAB  SSID  tool  and  it  is  the  only 

program  that  is  called  by  the  user.  MMLE  calls  the  supporting 

functions  to  provide  the  necessary  inputs  for  the  estimation 

algorithm.  The  primary  inputs  used  by  the  MMLE  macro  are; 

1)  uydata:  An  n-by-m  matrix  where  the  m  columns  are  the 
measured  time  history  inputs  of  the  vector  u  followed  by  the 
outputs  of  the  vector  y  and  n  is  the  number  of  time  steps. 
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2)  p2snam:  The  name  of  the  P2SS  function  that  converts  the 
parameter  vector  to  the  state-space  model. 

3)  pO:  The  initial  parameter  estimates  entered  during 

initialization . 

4)  rmsO:  A  vector  containing  the  standard  deviations  of  a- 
priori  parameter  estimates. 

5)  ggO:  The  initial  guess  of  the  innovations  co-variance 
matrix. 

6)  linesearch:  An  option  to  help  convergence  during  the 
Constrained-Newton  minimization  phase. 

7)  opt;  A  vector  that  specifies  the  niamber  of  iterations 
for  the  various  minimization  phases  and  the  convergence 
criteria  that  permits  termination  of  each  phase. 

For  this  study,  nasO  was  a  zero  vector  indicating  that 
initial  parameter  weighting  was  not  used.  The  initial 
innovations  co-variance  matrix,  ggO,  was  set  to  a  small,  non¬ 
zero  diagonal  matrix  to  minimize  the  effects  of  the  noise. 

There  are  four  variations  of  the  Newton  method  available 
for  the  minimization  of  the  weighted  error  sum.  The  opt 
vector  was  set  so  that  minimization  begins  with  a  Marquardt, 
gradient  algorithm.  The  gradient  algorithm  performs  well  when 
the  cost  function  is  very  steep,  but  performance  tends  to 
deteriorate  as  a  minimum  is  approached  and  the  cost  becomes 
more  flat  [Ref.  28;p.  SR-8].  When  the  user-specified  narq 
value  is  reached,  the  process  has  been  conditioned  well  enough 
for  a  Constrained-Newton  minimization  algorithm  to  begin. 
This  algorithm  is  most  accurate  near  the  minimum  of  a  cost 
function  when  the  cost  function  is  quadratic  [Ref.  9;p.  16]. 
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The  linesearch  option  is  set  "true"  to  help  convergence  in  the 
Constrained-Newton  phase.  Reference  28  discusses  each  of  the 
minimization  phases  in  more  detail.  For  this  study  the  value 
of  marq  was  set  at  0.02  and  the  convergence  criteria  was  set 
at  0.001.  All  parameters  are  identified  in  both  the  Marguardt 
and  the  Constrained-Newton  phases. 
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Rudder  Input  (deg)  Aileron  Input  (deg)  ^  I  Elevator  Input  (deg) 


BOM- 147  Elevator  Input 


r«  6  BQM-147  Longitudinal  Input 
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BOM-147  Ailenm  Input 
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Time  (secs) 


BOM-147  Rudder  Input 


Time  (secs) 


Flgura  7  BQM-147  Lateral-Directional  Input 


TABLE  3  FLIGHT  03-24  FLIGHT  CARD 


Pre-Flight 

and  Takeoff 

POl 

Engine  Off  Data 

P02 

Static  Measurement  (nose  up) 

P03 

Static  Measurement  (right  wing  up) 

P2A 

Takeoff  Data 

Stability  and  Control  Maneuvers  (100%) 

P2B1 

Pitch-Up  Doublet 

P2B2 

Pitch-Down  Doublet 

P2B3 

Roll-Right  Doublet 

P2B4 

Roll-Left  Doublet 

P2B5 

Yaw-Right  Doublet 

P2B6 

Yaw-Left  Doublet 

Trim  Flight 

Data 

P2E1 

20%  Throttle 

P2E2 

50%  Throttle 

P2E3 

100%  Throttle 
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TABLE  4  FLIGHT  03-24  MEASURED  SIGNALS 


Columns  1  through  4 


ALPHA 

BETA 

VELOCITY 

RPM 

Columns  5  through  8 

RUDDER 

AIL_LEFT 

ELEV_LEFT 

AIL_RIGHT 

Columns  9  through  12 

ELEV_RIGHT 

THROTTLE 

PITCH_RATE 

ROLLRATE 

Columns  13  through  16 

YAW_RATE 

AX 

AY 

AZ 

Columns  17  through  20 

AIL_CMD 

RUDDER_CMD 

ELEV_CMD 

ALTITUDE 

Columns  21  through  24 

SERVO  VOLT 


SERVO  AMP 


BLANK 


BLANK 


TABLE  5  BQM-147  PHYSICAL  CHARACTERISTICS 


Reference  Area  (S) 

21.24  ft^ 

Reference  Chord  (c-bar) 

4.458  ft 

Wingspan  (b-bar) 

8.167  ft 

Gross  Weight  (w) 

65  lb 

Moment  of  Inertia  (ly) 

3.443  slug-ft^ 

Moment  of  Inertia  (Ix) 

3.648  slug-ft^ 

Moment  of  Inertia  (Ixz) 

0 

Moment  of  Inertia  (Iz) 

5.965  slug-ft^ 

Sensor  Offset  (Xap) 

+2.5  ft 

Sensor  Offset  (Xbp) 

+2.5  ft 

Sensor  Offset  (Zan) 

0 

Sensor  Offset  (Zay) 

0 

Sensor  Offset  (Xan) 

0 

Sensor  Offset  (Xay) 

0 
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TABLE  6  INITIAL  PARAMETER  ESTIMATES,  LONGITUDINAL 


Parameter 

Value 

- 1 

Source  1 

Cl. 

2.4350 

Wind  Tunnel  [Ref  2:p  20] 

Cm. 

1.6900 

Wind  Tunnel  [Ref  2:p  20]  I 

Cm^ 

-0.1089 

DATCOM  [Ref  30]  I 

Clie 

-0.6670 

DATCOM  [Ref  30]  1 

C»ie 

-0.4060 

DATCOM  [Ref  30] 
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TABLE  7  INITIAL  PARAMETER  ESTIMATES,  LATERAL-DIRECTIONAL 


Parameter 

Value 

Source 

Cyfl 

-.2005  /rad 

Wind  Tunnel  [Ref  2;p  24] 

-.0325  /rad 

Smetana  [Ref  31 :p  108] 

CnB 

.0430  /rad 

Wind  Tunnel  [Ref  2:p  24] 

Clp 

-.3525  /rad 

Smetana  [Ref  31 :p  122] 

Cnp 

.0027  /rad 

Smetana  [Ref  31 :p  122] 

Cir 

.0729  /rad 

Smetana  [Ref  31:p  133] 

Cor 

-.0230  /rad 

Smetana  [Ref  31 ;p  135] 

.1377  /rad 

Smetana  [Ref  31 :p  124] 

-.0194  /rad 

Smetana  [Ref  31 :p  136] 

.1130  /rad 

Smetana  [Ref  31 :p  145] 

^l«r 

.0023  /rad 

Smetana  [Ref  31 :p  147] 

Cn«r 

-.0138  /rad 

Smetana  [Ref  31 :p  149] 

65 


VI .  pEst  APPLICATION 

For  many  aircraft  applications  the  linearized  equations  of 
motion  provide  a  reasonably  accurate  model  of  the  system  for 
small  disturbances  about  an  equilibrium  condition.  For  some 
applications,  however,  the  linear  model  is  insufficient  and  a 
more  accurate  non-linear  model  is  required.  Flight  testing  at 
the  Dryden  Flight  Research  Facility  highlighted  the 
limitations  of  the  linear  model  when  applied  to  maneuvers  at 
extreme  flight  conditions  and  with  some  unique  aircraft 
configurations  that  exhibited  non-linear  dynamic  behavior.  In 
1987  Murray  and  Maine  introduced  a  program  for  parameter 
estimation,  pEst,  in  non-linear  dynamic  systems.  A  full 
description  of  the  pEst  program  can  be  found  in  Reference  9. 

The  pEst  program  interfaces  directly  with  two  other 
programs.  The  GetData  program.  Reference  32,  is  a  FORTRAN 
utility  program  that  is  used  to  extract  selected  signals  from 
a  time-history  input  file  and  to  write  the  selected  signals  in 
a  user  specified  format  to  an  output  file.  The  second  program 
interfacing  with  pEst  is  XPLOT,  Reference  33.  The  XPLOT 
program  is  a  plotting  utility  that  is  used  to  plot  selected 
signals  from  the  measured  time-history  file  and  the  computed 
time-history  file.  The  following  sections  describe  how  the 
pEst  program  was  set  up  for  the  estimation  of  the  BQM-147 
stability  and  control  derivatives. 
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A.  MEASURED  TZME-HISTORY  DATA 

The  twenty-four  signals  that  comprise  the  recorded  time- 
history  data  from  the  flight  test.  Table  4,  were  re-formatted 
so  that  they  could  be  read  into  the  GetData  program.  The 
FORTRAN  utility  used  for  re-formatting  is  listed  in  Table  8. 


TABLE  8  DATA  CONVERSION  PROGRAM 


program  data 
dimension  dat(20) 

open(l,file='data01* , status= ' old ' ) 
k=0. 

read ( 1 , * ) 

20  read(l, 30,end=97)  (dat  (i),i=l,24) 

30  format(f6.0,2f7.0,f8.0,5f7.0,3x,f4.0,9f7.0, 

#  Ix,f7.0,f6.0,3f7.0) 
t=k/16. 

write(2,60)  t, (dat( j) , j*l,24) 

60  format(fl0.3,10x,3g20.14/5(4g20.14/) ,g20.14) 

k=k+l 
goto  20 

97  write(6,120) 

120  format (///'  ****  Normal  Termination  ***'///) 

99  stop 

end 


The  GetData  program  was  used  to  select  the  desired  signals 
from  the  available  time-history  data.  The  signals  not 
applicable  to  the  pEst  program  were  deleted.  Those  signals 
specifically  required  for  the  computation  of  the  estimates 
were  renamed  so  that  the  signal  names  on  file  matched  the 
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signal  names  expected  by  pEst.  One  new  signal  was  added  to 
provide  an  input  for  dynamic  pressure.  The  GetData  output 
file  was  saved  in  a  compressed  format  for  use  by  the  pEst 
program. 

B.  OPERATIONAL  STATUS 

The  pEst  program  uses  a  file  named  current  to  store  the 
operational  status  of  the  program.  The  current  file  is  used 
to  provide  initial  values  for  program  variables,  to  store  the 
status  of  the  program  for  recovery  and  to  store  the  results  at 
the  end  of  each  run.  The  primary  elements  of  the  current  file 
are  parameters,  states,  outputs  and  constants.  A  listing  of 
the  current  parameter  status  for  the  longitudinal  and  lateral- 
directional  applications  can  be  found  in  Tables  9  and  10. 

1 .  Parameters 

The  parameters  are  variables  that  can  be  estimated  by 
pEst.  The  five  attributes  of  each  parameter  are: 

1)  Current  value:  Typically  entered  by  the  user 

2)  Predicted  Value:  Computed  by  pEst 

3)  Estimation  Status:  Active  (T)  or  Inactive  (F) 

4 )  Cramer-Rao  bound 

5)  Change  in  value  from  previous  iteration 

The  current  values  of  active  parameters  were  the  same 
as  those  used  in  the  MMLE3  routine  multiplied  by  a  constant  to 
reflect  the  change  in  units  required  for  use  in  pEst.  The 
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active  parameters  were  those  that  are  listed  as  elements  of 
the  parameter  vectors  given  by  equations  (77)  and  (85) .  Only 
those  parameters  with  an  estimation  status  of  "true"  are 
estimated  by  the  program. 

The  parameter  list  also  allows  for  a  sensor  placement 
input.  In  the  MMLE3  routine  the  sensor  offset  was  corrected 
to  give  a  signal  that  simulated  a  measurement  at  the  eg.  In 
pEst,  the  effects  of  the  sensor  position  are  included  in  the 
observation  equations.  A  list  of  the  observation  equations 
can  be  found  on  page  20  of  Reference  9.  A  description  of  the 
available  parameters  can  be  found  in  Section  B  of  Reference  9. 

2 .  States 

The  status  of  the  state  variables  in  the  current  file 
is  used  to  determine  which  of  the  state  equations  are  to  be 
integrated.  The  state  equations  used  by  pEst  are  equations 
(35)  through  (43)  .  The  state  vectors  for  the  longitudinal 
model  and  the  lateral-directional  model  are  given  by  equations 
(75)  and  (82)  .  Since  roll  angle  was  not  measured,  the  state, 
phi,  was  set  inactive  for  the  lateral-directional  application. 

3 .  Outputs 

The  attributes  of  the  outputs  are  status  and 
weighting.  Default  weighting  was  used  in  the  initial  current 
file  in  order  to  avoid  biasing  the  solution  in  a  particular 
direction.  The  outputs  are  the  elements  of  the  output  vector, 
y,  of  the  form  given  by  equation  (48) . 
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4 .  Constants 


The  constants  are  variables  that  cannot  be  estimated 
by  pEst.  In  this  application,  constants  were  used  to  input 
the  physical  characteristic >  of  the  BQM-147.  The  constants 
and  their  respective  values  are  the  same  as  those  listed  in 
Table  5  for  the  MMLE3  application. 

C.  MINIMIZATION  ALGORITHMS 

The  pEst  program  has  several  minimization  algorithms 
available.  In  this  application  a  default  Newton-Raphson 
minimization  routine  was  used.  After  each  iteration  the 
active  parameters,  the  current  values  and  the  change  in 
current  valvis  from  the  previous  iteration  are  displayed.  The 
Cramer-Rao  bound  is  computed  with  the  Newton-Raphson  algorithm 
only.  The  Cramer-Rao  bounds  are  displayed  at  the  user's 
discretion. 

D.  COMPUTED  TIME-HISTORY 

For  each  output  whose  status  is  active  a  variable  time- 
history  is  computed.  The  'write'  command  is  used  to  write 
the  time  history  into  a  user  specified  output  file.  The 
signals  written  to  the  output  file  are  the  computed  states, 
the  computed  responses  and  the  residuals.  The  output  file  can 
be  read  by  the  XPLOT  program  and  user  selected  signals  can 
then  be  called  for  display. 
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TABLE  9  pEst  INITIAL  STATES;  LONGITUDINAL 


File  Name:  Current. p2bl  | 

Active  Parameters 

Value 

cNormO 

1.637171976 

cNorma 

0.046710000 

cNormde 

0.029500000 

cmO 

0.0 

cma 

-0.001900000 

cmq 

-0.677000000 

cmde 

-0.007100000 

alphaO 

0.0 

qO 

0.0 

ka 

1.0 

Constants 

Value 

mass 

2.0300 

ix 

3.6480 

iy 

3.4430 

iz 

5.9650 

ixy 

0.0 

iyz 

0.0 

ixz 

0.0 

area 

21.2400 

span 

8.1670 

chord 

4.4580 

States 

Status 

V 

F 

alpha 

T 

q 

T 

theta 

F 

Outputs 

Status 

V 

F 

alpha 

T 

q 

T 

theta 

F 

an 

T 
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TABLE  10  pEst  INITIAL  STATUS;  LAT-DIR 


File  Name:  Current. uav 

Active  Parameter 

Status 

cyO 

0.0 

cyb 

-0.00349900 

cydr 

0.00197210 

do 

0.0 

clb 

-0.00056719 

clp 

-0.35248000 

clr 

0.07290000 

clda 

0.00240300 

cldr 

0.00007888 

cnO 

0.0 

cnb 

0.00075044 

cnp 

-0.00270000 

cnr 

-0.02300000 

cnda 

0.00033857 

cndr 

-0.00055497 

kb 

1.0 

States 

S'i.at‘  f 

beta 

T 

P 

T 

r 

T 

phi 

T 

Outputs 

Status 

beta 

T 

P 

T 

r 

T 

phi 

T 

ay 

T 

Note:  Constants  are  the 

same  for  both 

1  Longitudinal  and  Lateral -Directional  | 

72 


VII.  RESULTS  AMD  DISCUSSION 

A  stated  purpose  of  this  study  was  the  practical 
application  of  parameter  estimation  for  estimating  stability 
and  control  derivatives  from  BQM-147  flight  test  data. 
Previous  information  from  a  wind-tunnel  study  was  made 
available;  however,  the  results  of  that  study  did  not  provide 
estimates  of  all  the  required  parameters,  and  changes  in  the 
production  configuration  would  be  expected  to  have  some 
influence  on  the  value  of  the  parameters.  Consequently, 
parameter  estimation  was  not  used  to  validate  wind-tunnel  or 
analytical  predictions  because  there  was  no  way  to  look  for 
consistency  between  the  data  sets,  nor  was  there  sufficient 
flight-test  data  for  a  valid  comparison.  The  complete 
aerodynamic  data-base  was  required  to  be  built  using  flight 
test  data  only.  The  wind-tunnel  and  analytical  predictions 
were  used  only  to  provide  starting  values  for  the  iterative 
algorithms  and  to  provide  a  reasonable  range  for  flight  test 
results. 

When  a  complete  data  base  must  be  built  using  flight  data 
only,  extensive  coverage  of  the  flight  envelope  is  needed. 
The  flight  test  data  made  available  for  this  study  comprised 
a  single  set  of  flight  maneuvers  consisting  of  pitch,  roll  and 
yaw  doublets.  For  the  longitudinal  case,  the  short-period 
mode  was  excited  by  an  elevator  doublet.  Table  3  shows  that 
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data  from  only  two  maneuvers,  a  pitch-up  doublet  (P2B1)  and  a 
pitch-down  doublet  (P2B2) ,  were  available  for  exciting  the 
short-period  mode.  The  lateral-directional  model  included  the 
roll,  dutch-roll  and  spiral  modes.  The  best  way  to  excite  all 
three  of  these  modes  is  to  have  both  aileron  and  rudder 
control  inputs.  Since  the  lateral-directional  maneuvers  were 
performed  under  similar  flight  conditions,  the  individual  roll 
(P2B3)  and  yaw  (P2B4)  doublets  were  paired  into  a  single 
maneuver.  The  measured  data  describing  the  maneuvers  and  the 
responses  are  found  in  Appendix  B. 

A.  MMLE3 

The  PC-based  MMLE3  routine  used  in  this  application  relies 
on  the  locally  linearized  aerodynamic  model  presented  in 
Chapter  V.  The  results  of  the  wind-tunnel  study  [Ref.  2]  show 
that  the  behavior  of  the  vehicle  is  approximately  linear  for 
angles  of  attack  up  to  about  15  degrees.  The  linearized  model 
is  valid  only  when  the  assumptions  made  in  the  linearization 
process  are  valid.  Thus,  the  model  is  ideally  suited  for 
small  angle  perturbation  maneuvers  and  it  can  be  expected  that 
as  the  magnitude  of  the  maneuver  increases,  the  locally 
linearized  model  loses  its  validity  and  the  estimates  from  the 
model  are  a  less  accurate  representation  of  the  aircraft 
dynamics. 

The  process  noise  option  was  used  in  tae  MMLE3  routine. 
The  risk  in  using  the  process  noise  option  is  that  the  Kalman 
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gains  tend  to  become  high,  giving  the  high-frequency  modeling 
errors  more  opportunity  to  influence  parameter  values. 

1.  Model  Validation 

In  Chapter  III,  the  changes  that  were  necessary  to 
Graham's  routine  [Ref.  6]  to  reduce  the  order  of  the  model 
were  discussed.  In  order  to  validate  the  changes,  the 
modified  model  was  run  using  the  same  F-14A  data  that  were 
presented  in  Reference  6.  The  results  of  the  run  shown  in 
Figures  8  through  11  compare  well  with  Graham's  results. 

2 .  Longitudinal  Model 

The  computations  in  the  MMLE3  routine  were  limited  by 
the  working  precision  of  MATLAB.  Prior  experience  showed  that 
the  PC-based  routine  was  highly  sensitive  to  the  initial  input 
derivatives.  After  several  adjustments  to  the  initial 
parameter  vector,  the  run  using  the  P2B1  maneuver  converged 
successfully.  The  run  using  the  P2B2  maneuver  did  not 
converge  because  of  a  singularity  due  to  working  precision. 
The  results  of  the  P2B1  maneuver  are  shown  in  Figures  12 
through  15.  The  MATLAB  output  listing  showing  the  input  and 
determined  derivatives,  per  radian,  is  found  in  Table  11. 

The  estimated  response  for  a  appears  to  lag  the  measured 
response  by  a  time  step  at  the  beginning  of  the  maneuver  and 
is  unable  to  keep  up  with  the  amplitude  of  the  measured 
response  in  the  pitch-down  portion  of  the  maneuver.  This 
behavior  could  have  been  predicted  since  pitch  angle,  which 
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was  not  measured,  is  an  important  term  in  the  a  state  equation 
when  the  maneuver  is  not  restricted  to  small  perturbations 
from  a  stable  flight  condition.  The  estimated  q  response 
generally  correlates  well  with  the  measured  q  response  though 
it  does  not  pick  up  the  high  rates  at  the  extreme  ends  of  the 
maneuver.  The  differences  in  q  may  be  due  to  the  ignored 
terms  in  the  state  equation  in  the  linearization  process.  The 
estimated  normal  acceleration  response  appears  to  correlate 
well. 

3 .  Lateral-Directional  Model 

Inspection  of  the  data  in  Appendix  B  indicates  that 
the  data  were  exceptionally  noisy  in  the  absence  of  any 
control  input.  Noise  in  the  general  frequency  range  of  the 
system  response  was  a  problem.  Vfhen  the  P2B3  and  P2B5 
maneuvers  were  paired  for  the  lateral-directional  run,  the 
process  noise  was  too  large  and  the  Gauss-Newton  iteration  was 
terminated.  Consequently  there  were  no  results  for  the 
lateral-directional  model  in  the  MMLE3  routine. 
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Figure  9  Model  Validation:  F-14  AOA  Response 
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Figur*  12  BQM-147  Elevator  Input  (P2B1) 
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Figure  15  BQM-147  Normal  Acceleration  Response  (P2B1) 
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TABLE  11  MATLAB  RESULTS  FOR  P2B1 


NFS  PARAMETER  IDENTIFICATION  MACRO  FILE 
FOR  ACTUAL  FLIGHT  TESTS  USING  SIMPLFIEC  SHORT  PERIOD 
LONGITUDINAL  STABILITY  AND  CONTROL  DERIVATIVES 


RESULTS 

pid 

p(pid) 

pref 

cramer 

2 f cramer 

insens 

gdop 

1.0 

2.7391 

2.4350 

0.1084 

0.3085 

0.0780 

1.3902 

2.0 

-0.0887 

-0.1089 

0.0028 

0.0079 

0.0018 

1.5088 

3.0 

-0.1078 

-0.6770 

0.0414 

0.1178 

0.0159 

2.6088 

4.0 

0.5015 

1.6900 

0.0763 

0.2171 

0.0666 

1.1457 

5.0 

-0.0683 

-0.4060 

0.0042 

0.0120 

0.0017 

2.4701 

MMLE 

STABILITY 

&  CONTROL 

DERIVATIVES 

CLA 

CMA 

CMQ 

CLDE 

CMDE 

2.7391 

-0.0887 

-0.1078 

0.5015 

-0.0683 

INITIAL 

INPUT  DERIVATIVES 

2.4350 

-0.1089 

-0.6770 

1.6900 

-0.4060 
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B.  pBst 

1.  Model  Validation 

A  secondary  purpose  of  this  study  was  to  integrate  the 
pEst  parameter  estimation  routine  into  the  ongoing  research 
program  at  the  Naval  Postgraduate  School,  Department  of 
Aeronautics  and  Astronautics.  Successful  test  cases  for  both 
the  longitudinal  and  lateral-directional  models  were  run  to 
validate  the  model.  The  pEst  parameter  estimation  routine 
should  provide  greater  capability  to  the  Department  because  it 
can  be  used  to  model  the  dynamic  behavior  of  aircraft  that 
cannot  be  appropriately  modeled  using  the  shorter  PC-based 
routine. 

2 .  Longitudinal  Model 

The  longitudinal  case  was  validated  using  a  T-37  time 
history  file  received  from  NASA  Dryden.  The  results  of  the 
validation  are  presented  in  Figures  16  through  19.  The 
initial  set  up  for  the  BQM-147  longitudinal  model  was 
presented  in  Table  9  in  Chapter  VI.  The  same  current  file  was 
used  for  both  the  P2B1  and  P2B2  maneuvers. 

a.  P2B1 

The  results  of  the  P2B1  maneuver  are  shown  in 
Figures  20  through  23  and  the  pEst  output  listing  is  shown  in 
Table  12.  The  estimated  response  for  a  correlates  well  except 
at  the  extreme  edge  of  the  pitch-down  part  of  the  maneuver. 
The  measured  data  shows  a  sinusoidal  response  at  the  beginning 
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of  the  window  prior  to  any  control  input.  This  may  be  an 
accentuation  of  the  process  noise  due  to  the  low  sample  rate. 
On  the  initial  run,  the  Cramer-Rao  bound  for  was 

unacceptably  high.  Rather  than  accept  the  high  bound,  weight 
estimates  based  on  discussions  with  NASA  Dryden  personnel, 
were  applied  to  the  outputs  and  subsequent  runs  showed 
improved  results.  The  adjustment  gave  a  better  estimated 
response  fit  for  a  and  a  more  reasonable  Cramer-Rao  bound,  but 
the  tradeoff  was  an  increase  in  the  pitch-damping  coefficient 
which  is  evident  in  the  q  response.  The  response  presented 
for  both  a  and  q  represents  the  best  possible  fit  achieved 
from  multiple  combinations  of  output  weighting.  The 

estimated  response  for  normal  acceleration  correlates  well, 
though  it  seems  to  lead  the  measured  response  by  a  time  step. 
No  time  skewing  was  attempted. 
b.  P2B2 

The  results  of  the  P2B2  maneuver  are  shown  in 
Figures  24  through  27,  and  the  pEst  output  listing  is  shown  in 
Table  13.  The  final  output  weighting  values  used  in  the  P2B1 
case  were  used  for  the  P2B2  case  because  they  appeared  to  give 
the  best  combination  of  results.  The  estimated  a  response 
does  not  keep  up  with  the  magnitude  of  the  measured  response 
in  the  pitch-down  part  of  the  maneuver  and  tends  to  dampen 
more  quickly  at  the  end  of  the  maneuver.  Again,  this  may  be 
due  to  ignoring  the  ©  term  in  the  a  state  equation.  The 
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measured  data  show  the  same  sinusoidal  response  in  the  absence 
of  control  input  at  the  beginning  of  the  time  window  and  again 
at  the  end.  The  estimated  q  response  correlates  well.  The 
damping  on  q  tended  to  increase  as  weight  on  a  was  increased 
indicating  a  corellation  between  the  a  term  and  the  q  term. 
Normal  acceleration  fits  well  though  it  does  not  pick  up  the 
peaks  at  the  edges  of  the  maneuver. 

3 .  Lateral-Directional  Model 

The  lateral-directional  case  was  validated  using  a  PA- 
30  time  history  file  received  from  NASA  Dryden.  The  results 
of  the  validation  are  shown  in  Figures  28  through  33.  The 
results  of  the  combined  P2B3-P2B5  lateral-directional  maneuver 
are  shown  in  Figures  34  through  39.  The  pEst  output  listing 
is  presented  in  Table  14.  The  6  response  correlates  well  at 
the  onset  of  the  aileron  input  and  follows  the  measured 
response  through  the  rudder  doublet.  The  r  response  fits  well 
at  the  onset  of  the  aileron  input;  however,  the  computed 
response  appears  to  dampen  more  quickly  than  the  measured 
response  at  the  end  of  the  rudder  doublet.  The  p  response 
does  not  follow  the  rate  chrough  the  roll-left  portion  of  the 
aileron  doublet.  The  estimated  response  fits  well  during  the 
rudder  input  but  dampens  more  rapidly  than  the  measured 
response.  The  roll  angle  response  does  not  correlate  well 
with  the  approximated  angle  data  obtained  through  numerical 
integration  of  the  roll  rate.  The  lateral  acceleration 
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response  seems  to  match  in  phase  but  not  in  amplitude.  The 
measured  response  shows  better  damping  at  the  completion  of 
the  rudder  doublet. 

4 .  Summary 

A  summary  of  the  longitudinal  parameters  is 
presented  in  Table  15.  The  columns  labeled  P2B1  and  P2B2 
contain  the  results  of  the  pEst  parameter  estimation  runs. 
The  Cramer-Rao  bounds  represent  an  approximation  to  the 
standard-deviation  and  are  consistently  one  order  of  magnitude 
less  than  the  value  of  the  parameter.  The  large  perturbation 
input  maneuver  and  subsequent  responses  should  lead  one  to 
suspect  that  the  linear,  PC  based  MMLE3  routine  would  be  an 
insufficient  model  for  this  case  study.  The  plots  shown  in 
Figures  12  through  15  support  this  suspicion.  The  pEst 
results  provide  a  better  fit  between  the  measured  and 
estimated  data.  The  confidence  level  is  higher  for  the 
parameters  Cm,  and  Cm,,,  which  have  lower  relative  Cramer- 
Rao  bounds.  The  values  for  C^,  differ  by  a  factor  of  2  and 
the  sign  on  Cm^  was  expected  to  be  negative,  indicating  a 
restoring  moment  for  positive  pitch  rate.  The  positive  Cm^ 
may  indicate  that  there  is  a  strong  correlation  between  the  a 
term  and  the  q  term.  More  data  points  are  required  to  improve 
the  confidence  level  for  these  two  parameters. 

The  lateral-directional  case  is  summarized  in  Table  16. 
Running  the  lateral-directional  case  in  the  non-linear  PC 
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based  routine  was  unsuccessful  because  the  process  noise  was 
too  high.  Only  the  pEst  results  are  shown.  The  B  derivatives 
and  the  p  derivatives  have  reasonable  Cramer-Rao  bounds.  The 
variance  in  the  r  derivatives  is  unacceptably  high  and 
indicates  a  low  confidence  level  in  the  value  of  the 
parameters.  A  sign  change  in  some  of  the  control  derivatives 
indicates  a  strong  correlation  with  damping  terms,  a  noted 
problem  with  parameter  estimation  methods.  More  flight  data 
are  required  for  any  further  analysis. 

Three  options  are  available  to  correct  for  the  high 
Cramer-Rao  bound. 

1)  Fix  the  value  of  the  parameter 

2)  Weight  the  estimate 

3)  Use  a  better  input  stimuli 

Without  a-priori  knowledge  of  the  parameters,  fixing  the 
value  of  the  parameter  was  not  a  practical  option.  The 
weighting  option  did  not  help  to  improve  the  overall  results 
because  of  the  tradeoffs  in  the  values  of  other  parameters. 
A  more  reasonable  approach  to  output  weighting  could  have  been 
made  only  if  there  were  data  other  than  the  single  flight  test 
supporting  the  aerodynamic  data  base.  With  the  limited 
maneuvers,  better  input  stimuli  were  not  available.  Without 
any  of  the  options  available,  the  only  choice  is  to  accept  the 
high  Cramer-Rao  bound  with  the  understanding  that  the 
confidence  level  in  the  value  of  the  parameter  is  very  low. 
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There  were  three  limitations  in  the  data  that  affected  the 
outcome  of  the  lateral-directional  results.  Signal  coverage 
is  very  important  and  unless  the  flight  test  maneuvers  are 
restricted  to  small  perturbations  about  a  stable  flight 
condition,  the  body  angles  should  be  included  in  the  signal 
set.  Noise  in  the  same  general  freguency  range  as  the 
response  is  a  problem.  The  flight  data  were  sampled  at  16  Hz 
and  nothing  is  known  about  the  sensor  accuracies,  filtering, 
digitizing  and  time-tagging  the  data.  The  low  sample  rate  may 
have  distorted  the  response  signal  by  accentuating  the  process 
noise.  To  excite  the  lateral-directional  modes,  rudder  and 
aileron  doublets  should  be  included  in  a  single  input 
maneuver,  separated  by  sufficient  time  to  allow  for  the  first 
response  to  dampen.  In  Flight  03-24,  the  maneuvers  were 
independent  and  combining  the  data  sets  together  may  have 
induced  additional  error  in  the  response  estimates. 
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re  18  pEst  Validation:  q  Response 


pE<it  Validation  Normal  Acceleration 


Figure  19  pEst  Validation:  Normal  Acceleration 
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Figure  21  BQM-147  a  Response  (P2F1) 
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Figure  23  BQM-147  Normal  Acceleration  (P2B1) 
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TABLE  12  pEst  OUTPUT  LISTING 


Iteration 

Lev-Harq 

17 

used  1  trials. 

V  •  0 . 

len  * 

3.150 

cost  ■  37.65 

name 

value  on? 

delta 

name 

value 

on? 

delta 

cNormO 

1.173  T  - 

.4573B-02 

o&q 

0.3234 

T 

0.2962E-03 

cNorma 

0.8794E-01  T  0 

.r699E-03 

cmde 

- .1466E- 

-02  T 

0.4741E-05 

cNormde 

0.3389E-01  T  - 

.3382E-03 

alphaO 

-.5625 

T 

-.1678E-02 

cmO 

-.3865E-01  T  - 

.1539E-05 

ka 

1.663 

T 

0.8075E-02 

cma 

-.2951E-02  T  - 

.1587E-04 

total  cost  >  37.65 

cost  per  response  : 

alpha  q  an 

35.5  69.5  7.88 

Iteration  17  used  12  integrations 
***  iteration  converged 
estimate  used  209  integrations 


pEst:  par 
name 

♦cr 

value  on? 

cr  bound 

name 

value 

on? 

cr  bound 

cMormO 

1.173 

T 

0.4632E-01 

cmq 

0.3234 

T 

0.2749E-01 

cNorma 

0.8794E-01 

T 

0.3379E-02 

cmde 

-.1466E 

■02  T 

0.413eE-04 

cNormda 

0.23e9E-01 

T 

0.2481E-02 

alphaO 

-  .5625 

T 

0.1954 

cmO 

-.3865E-01 

T 

0.1263E-02 

ka 

1.663 

T 

0.6568E-01 

cma 

pEst; 

- .2951E-02 

T 

0.1338E-03 
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Figure  25  BQM-147  a  Response  (P2B2) 
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Figure  27  BQM-147  Normal  Acceleration  (P2B2) 
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TABLE  13  pEst  OUTPUT  LISTING  (P2B2) 


Iteration  14 


Lev-Marc  a 

sed  1  trials 

.  V  «  0. 

l«n  « 

1.637  cost  .  103.0 

name 

value  on? 

delta 

name 

value  on? 

delta 

cNormO 

0.6941  T 

-.1956E-02 

cma 

-.2266E-02  T 

-  .3263E-05 

cNorma 

0.4120E-01  T 

-.1341E-03 

C9q 

0.1355  T 

0.3866E-02 

cNormde 

0.2149E-01  T 

-.5021E-04 

cmde 

-.1600E-02  T 

0.1210E-04 

emO 

-.3293E-01  T 

0.2963E-05 

lea 

1.573  T 

-  .4011E-04 

total  cost 

.  103.0 

cost  per  response  : 

alpha  q  an 

6«.9  223.  19.5 

Iteration  14  used  11  integrations 
•»*  iteration  converged 
estimate  used  157  integrations 

pEst:  write  out2.p2b2 

writing  cemputed  time  history  to  file  *out2.p2b2* 
pEst;  par  +cr 


name 

value  on? 

cr  bound 

naaa 

value  on? 

cr  bound 

cNormO 

0.6941 

T 

0.2587E-01 

CM 

-.2266E-02 

T 

0.8947E-04 

cNorma 

0.4120E-01 

T 

0.1786E-02 

oaq 

0.1355 

T 

0.1741E-01 

cNormde 

0.2149E-01 

T 

0.1669E-02 

ciada 

-  .1600E-02 

T 

0.4203E-04 

emO 

- .3293E-01 

T 

0. 1119E-02 

ka 

1  .573 

T 

0 .3094E-01 
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Roll  Rate  (deg/sec)  C  Yaw  Rate  (deg^) 


pEst  Validation  Yaw  Rate 


re  30  pEst  Validation;  r  Response 


pEst  Validation  Roll  Rate 


Figure  31  pEst  Validation:  p  Response 
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Figure  34  BQM-147  Lateral-Directional  Input 
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Figure  35  BQM-147  B  Response 
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Roll  Rate  (de0ec)  ^  |  Yaw  Rate  (deg/Mc) 


Flgur*  37  BQM-147  p  Response 
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Figure  38  BQM-147  (p  Response 


101 


TABLB  14  pEst  OUTPUT;  LATERAL-DIRECTIONAL 


Lev-Narq 

used  1  trials,  v  >  0. 

len  » 

3.436  cost  «  1003. 

name 

value  on? 

delta 

name 

value  on7  delta 

cyO 

0.4222E-02 

T 

-.5625E-06 

cldr 

0.38S6E-04 

T  -.2546E-06 

cyb 

-.8251E-02 

T 

-  .n25E-04 

cnO 

0.3708E-03 

T  -.2122E-05 

cydr 

-.1882E-02 

T 

0.2043E-OS 

cnb 

0.4335E-03 

T  0.3397E-06 

elO 

0.1080E-02 

T 

0.3836E-0S 

cnp 

-.1212E-01 

T  0.1558E-03 

clb 

-.9698E-04 

T 

0.2407E-Q7 

cnr 

-.7893E-03 

T  0.1786E-03 

clp 

-.5860E-01 

T 

-.2733E-03 

cnda 

-.7178E-04 

T  0.3820E-06 

clr 

-.3279E-02 

T 

-.1693E-03 

crdr 

-.2109E-03 

T  0.2517E-06 

clda 

-.S099E-03 

T 

-  .1818E-0S 

kb 

1 .067 

T  0.2492E-02 

CoCal  coac  •  1003. 

coaC  per  response  i 

beta  p  r  phi 

51.2  3.714E»03  322.  817 

Iteration  8  used  19  integrations 
iteration  converged 
estimate  used  155  integrations 

pEst:  par  4cr 


name 

value  on7 

cr  bound 

name 

value  on? 

cr  bound 

cyO 

0.4222E-02 

T 

0.272SE-02 

cldr 

0.3856E-04 

T 

0  8644E-05 

cyb 

-.82S1E-02 

T 

0.1162E-02 

cnO 

0.3708E-03 

T 

0.S76SE-04 

cydr 

-.1882E-02 

T 

0.6809E-03 

cnb 

0.4335E-03 

T 

0.1105E-04 

CIO 

0.1080E-02 

T 

0.4500E-04 

cnp 

-  .1212E-01 

T 

0.3082E'02 

Clb 

-  .9698E-04 

T 

0.9580E-05 

cnr 

-  .7893E-03 

T 

0.6030E-02 

clp 

-  .5860E-01 

T 

0.2814E-02 

cnda 

-  .7178E-04 

T 

0. 180SE-04 

clr 

-.3279E-02 

T 

0.2653E-02 

cndr 

-  .2109E-03 

T 

0.1095E-04 

clda 

-.6099E-03 

T 

0.2033E-04 

kb 

1.067 

T 

0.9189E-01 

pEst:  write  outlat7 


»y 

109. 
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TABLE  15  LONGITUDINAL  SUMMARY 


Parameter 

Initial 

MMLE3 

P2B1 

P2B2 

Cta 

2.4350 

2.7391 

5.0367 

2.3600 

CLfe 

1.6900 

0.5015 

1.3689 

1.2314 

Cm. 

-0.1089 

-0.0887 

-0.1691 

-0.1298 

Cmq 

-0.6670 

-0.1078 

0.4997 

0.1355 

-0.4060 

-0.0683 

-0.0840 

-0.0917  1 

Cramer-Rao  Bounds 

Cl. 

0.1084 

0.1936 

0.1023 

CL^. 

0.0763 

0.1422 

0.1071 

Cm. 

0.0028 

0.0077 

0.0051 

Cm^ 

0.0414 

0.0275 

0.0174 

CTOae 

0.0042 

0.0024 

0.0024 
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TABLE  16  LATERAL  DIRECTIONAL  SUMMARY 


Parameter 

Initial 

pEst 

Cramer-Rao  | 

CyB 

-0.2005 

-0.4728 

0.0665 

^16 

-.0325 

-0.0056 

0.0005 

Cng 

0.0430 

0.0248 

0.0006 

Clp 

-0.3525 

-0.0586 

0.0028 

Cnp 

0.0027 

-0.0121 

0.0031 

cip 

0.0729 

0.0033 

0.0026 

Cn, 

-0.0230 

-0.0008 

0.0060 

Clia 

0.1377 

-0.0349 

0.0012 

-0.0194 

-0.0041 

0.0010 

0.1130 

-0.1078 

0.0390 

Cl«r 

0.0023 

0.0022 

C.0005 

-0.0138 

-0.0121 

0.0006  1 
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VIII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Based  on  the  results  presented  in  Chapter  7,  the  following 
conclusions  are  made. 

1)  The  PC-based  parameter  estimation  routine  is  ideally 
suited  for  small  maneuvers  on  stable  platforms.  The 
magnitude  of  the  maneuvers  in  Flight  03-24  exceeded  the 
range  of  validity  for  the  linear  model. 

2)  The  pEst  routine  was  successfully  integrated  into  the 
NPS  Department  of  Aeronautics  and  Astronautics  flight 
research  program. 

3)  Noise  in  the  same  general  frequency  range  as  the  system 
response  constitutes  a  major  problem. 

4)  Body  angles  should  be  included  in  the  signal  set  unless 
the  maneuvers  are  restricted  to  small  perturbations  about  a 
stable  flight  condition. 

5)  Improving  the  confidence  level  in  the  values  of  the 
parameters  presented  in  Tables  15  and  16  will  require 
modifications  to  the  flight  test  program. 


B.  RECOMMENDATIONS 

Based  on  the  conclusions  presented  above  and  the  problems 
encountered  during  the  course  of  this  study,  the  following 
recommendations  are  made. 

For  the  Marine  Corps  UAV  Project  Office,  Intelligence 
Systems  (C2IU) ,  Marine  Corps  Research,  Development  and 
Acquisition  Command,  Quantico,  VA.  : 
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1)  Evaluate  the  current  status  of  the  BQM-147  Program  and 
determine  if  additional  flight  testing  to  provide  higher 
confidence  values  of  stability  and  control  parameters  is 
warrant  *»d. 

2)  If  higher  confidence  values  are  desired,  formulate  a 
clear  set  of  objectives  and  modify  the  flight  test  program 
to  include  a  series  of  multiple  maneuvers  designed 
specifically  to  meet  those  objectives.  The  maneuvers  should 
include  repeated  data  points,  and  roll  and  pitch  angles 
should  be  measured. 


For  the  NPS,  Department  of  Aeronautics  and  Astronautics: 

1)  Update  the  Department  PC  Lab  to  include  MATLAB  version 
3.5j  or  higher.  As  a  part  of  the  update,  add  the  State- 
Space  Identification  Tool  to  the  current  inventory  of  MATLAB 
toolboxes.  The  update  will  allow  instructional  use  of  the 
PC-based  routine  in  the  flight  dynamics  and  systems  courses 
and  practical  application  of  a  desk-top  routine  at  the  NPS 
UAV  Flight  Research  Lab. 

2)  Incorporate  pEst  as  an  instructional  tool  in  the  Flight 
Dynamics  courses  at  NPS. 

3)  Investigate  modification  of  the  user-modifiable 
subroutines  in  pEst  to  include  helicopter  dynamics  and 
missile  dynamics. 

4)  Investigate  the  feasibility  of  offering  a  course  in  non¬ 
linear  flight  dynamics  or  non-linear  control  of  aerospace 
vehicles. 

5)  Suggest  as  a  potential  project  for  the  Avionics  System 
Design  or  Aircraft  Design  courses,  the  design  of  an 
instrumentation  package  for  a  UAV  Flight  Test  Program  that 
would  respond  to  a  real-world  DoD  requirement. 
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Flgur*  A-1  EXDRONE  Baseline  Configuration 
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AOA  (deg)  •©  Elevator  Input  (deg) 


APPENDIX  B  FLIGHT  03-24  TIME  HISTORY  DATA 


A.  PITCH-UP  DOUBLET  (P2B1) 


Figure  B-2  Angle  of  Attack  Response 


113 


AOA  (deg)  I  g  Elevator  Input  (deg) 


PITCH>DOWN  DOUBLET 


re  B-5  Elevator  Control  Input 


Figure  B-6  Angle  of  Attack  Response 
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C.  COMBINED  ROLL  AND  YAW  DOUBLET  (P2B3-P2B5) 


Figure  B-10  Sideslip  Response 
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RoU  Rate  (deg/tec)  «g  Yaw  Rate  (deg/sec) 


Lateral  Accelention  (g) 


Figure  B-13  Lateral  Acceleration 
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